99 


IEP€  99-217 

Computation  of  Current  to  a  Moving  Bare  Tether 

T.  Onishi  and  M,  Martinez- Sanchez, 

MIT,  Cambridge,  MA  USA 

D.L.  Cooke, 

Ait  Force  Research  Laboratory, 

Space  Vehicle  Directorate,  Hanscom  AFB, 

MA,  USA 


26th  International  ICIectriej  Propulsion  Conference 
October  17-21, 1999,  Kitakyuslni,  Japan 


REPORT  DOCUMENTATION  PAGE 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  inc 
gathering  and  maintaining  the  data  needed,  and  completing  and  r^iewing  the  collection  of  information, 
collection  of  information,  including  suggestions  for  reducing  this  t-wden,  to  Washington  Headquarters  Se 
Davis  Highway,  Suite  1  204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  P 


Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4 

1.  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE 


REPORT  DATE 
18  Dec  00 


Budget,  P 

I  3.  RE 


afrl-sr-bl-tr-00- 


Final  Report:  15  Jul  95  To  14  Jul  99 


5.  FUNDING  NUMBERS 


3-D  Model  Developments  for  Plasma  Plumes 


F49620-95- 1-0444 


6.  AUTHOR(S) 

Prof.  D.E.  Hastings 
Prof.  M.  Martinez  Sanchez 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Massachusetts  Institute  Of  Technology 
77  Massachusetts  Avenue 
Room  El 9-7 19 

Cambridge,  Massachusetts  02139 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING/MONITORING  AGENCY  NAME{S)  AND  ADDRESS(ES) 

AFOSR 

801  N.  Randolph  St. 

Arlington  VA  22203 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  Public  Release:  Distribution  is  unlimited 


13.  ABSTRACT  (Maximum  200  words ) 

Two  different  studies  were  supported  under  this  AASERT  Grant,  One,  performed  by  Master's  Candidate,  Gregory  Yashko, 
examined  the  technical  feasibility  and  economic  impact  of  clusters  of  micro-satellites  propelled  by  micro-ion  engines.  It  was 
concluded  that  sensor  distribution  among  several  satellites  is  appropriate  when  savings  from  the  increased  reliability  outweigh 
the  increased  deployment  costs,  and  that  large  constellations  of  g-satellites  could  feasibly  perform  missions  currently  being 
carried  out  by  large  satellites.  It  was  also  found  that  the  performance  of  micro-fabricated  linear  ion  thrusters  ~s  poor  because 
of  insufficient  primary  electron  confinement,  so  other  micro-propulsion  concepts  appear  to  be  necessary  for  implementation 
of  g-sat  clusters.  The  other  study  was  carried  out  by  Master's  (later  Ph.D.)  candidate,  Tatsuo  Onishi,  and  it  concerned  the 
numerical  modeling  of  electron  capture  by  an  ionospheric  orbiting  bare  tether.  A  PIC  method  was  developed  for  the 
simplified  case  in  which  magnetic  field  and  relative  plasma  motion  are  ignored,  and  good  agreement  was  found  with  classical 
results  of  Probe  Theory  for  this  case.  These  results  opened  the  way  for  simulation  of  realistic  tether  configurations,  a  work 
which  is  now  continuing  under  separate  funding. 


14.  SUBJECT  TERMS 


3-D  Model  Developments  For  Plasma  Plumes 


15.  NUMBER  OF  PAGES 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRAC  I 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


UNCLASSIFIED 


UNCLASSIFIED 


UNCLASSIFIED 


Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


Computation  of  Current  to  a  Moving  Bare  Tether 

Tatsuo  Onishi*and  Manuel  Martinez-Sanchez^ 

Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts,  USA 

David  L.  Cooke 

Air  Force  Research  Laboratory 
Space  Vehicle  Directorate ,  Hans  com  AFB,  MA,  USA 


Abstract 

An  electrodynamic  bare  tether  has  been  consid¬ 
ered  as  an  alternative  method  of  propulsion  with¬ 
out  expenditure  of  propellant. The  object  of  the  work 
reported  here  is  the  development  of  a  numerical 
method,  Particle-In- Cell  method,  for  the  calculation 
of  electron  current  to  a  positive  bare  tether  moving 
at  orbital  velocity  in  the  ionosphere,  i.e.  in  a  flow¬ 
ing  magnetized  plasma  under  Maxwellian  collision¬ 
less  conditions.  The  code  uses  the  quasi-neutrality 
condition  to  solve  for  the  local  potential  at  points  in 
the  plasma  which  coincide  with  the  computational 
outside  boundary.  Given  the  boundary  conditions, 
Poisson  equation  is  solved  in  such  a  way  that  the 
presheath  region  can  be  captured  in  the  computa¬ 
tion.  Electrons  are  assumed  to  have  a  Maxwellian 
distribution  at  the  boundary  due  to  their  high  mo¬ 
bility,  whereas  ions  are  assumed  so  only  in  the  far 
upstream  region  and  are  also  assumed  to  only  decel¬ 
erate  one  dimensionally  due  to  their  large  mass.  The 
results  indicate  a  stable  convergence,  and  clearly 
represent  a  presheath  region.  Collected  currents 
turn  out  to  be  more  than  that  predicted  by  the  Or¬ 
bital  Motion  Limit  (OML)  theory. 

Introduction 

In  tethered  satellite  technology,  it  is  important  to  es¬ 
timate  how  many  electrons  a  spacecraft  can  collect 
from  its  ambient  plasma  by  its  positively  charged 
electrodynamic  tether.  The  analysis  is,  however, 
quite  complicated  because  of  the  small  but  signifi¬ 
cant  Geomagnetic  field  and  the  spacecraft’s  relative 
motion  to  both  ions  and  electrons.  One  of  the  ap¬ 
proaches  to  this  solution  is  the  numerical  method.  In 
the  numerical  analysis  of  space  plasma,  one  of  the 
most  reliable  methods  has  been  the  Particle-In- Cell 
(PIC)  method.  In  this  paper  we  develop  a  PIC  code 
for  a  two  dimensional  collisionless  plasma  under  the 
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effects  of  magnetic  field  and  the  spacecraft’s  relative 
motion. 

This  paper  is  the  continuation  of  the  authors’  work 
on  a  quiescent  unmagnetized  plasma  [1].  In  [1],  the 
quasi-neutral  condition  was  used  to  obtain  the  out¬ 
side  boundary  condition  for  the  Poisson  solver.  The 
quasi-neutral  codition  requires  local  densities  of  elec¬ 
trons  and  ions  at  points  on  the  outside  boundary.  In 
computation,  each  density  consists  of  two  parts;  one 
due  to  the  incoming  particles  and  the  other  due  to 
the  outgoing  particles.  Based  on  same  assumptions, 
the  density  of  the  incoming  particles  can  be  analyzed 
and  expressed  as  a  function  of  the  local  potential. 
The  density  of  the  outgoing  particles  is  computed 
by  using  the  same  technique  in  [1],  namely,  by  esti¬ 
mation  from  the  numerical  results  for  the  previous 
time  step. 

Since  some  electrons  are  absorbed  on  the  surface 
of  the  tether,  there  is  also  a  region  called  “presheath” 
outside  the  sheath,  where  the  quasi-neutrality  pre¬ 
vails  but  the  electric  potential  is  not  equal  to  that 
of  the  ambient  plasma  at  infinity.  Thus,  in  order 
to  reproduce  the  presheath  region  in  the  computa¬ 
tion,  the  quasi-neutral  condition  is  physically  very 
reasonable. 

In  the  wake  region  behind  the  tether,  the  Debye 
lengths  of  electrons  become  very  large  due  to  the  low 
electron  density,  which  in  turn  is  due  to  the  fact  that 
ions  are  deflected  by  the  tether  potential  and  thus 
do  not  reside  in  that  region.  The  large  Debye  length 
may  exclude  the  quasi-neutral  condition  where  the 
computational  grid  is  smaller  than  the  local  Debye 
length.  In  our  computation,  we  apply  the  quasi¬ 
neutrality  where  we  can  obtain  the  potential  from 
it,  otherwise  we  use  the  zero-gradient  condition. 

Computation 

The  major  difficulty  of  a  PIC  method  applied  to 
an  infinitely  large  plasma  appears  in  the  spec¬ 
ification  of  the  computational  outside  boundary 
condition,  namely  the  velocity  distribution  func- 


Magnetic  field 

0.3  Gauss 

Ion  mass  (0+) 

2.67  x  nr27  kg 

Electron  temperature 

0.1  eV 

Ion  temperature 

0.1  eV 

Electron  thermal  velocity 

212  km/sec 

Ion  thermal  velocity 

1  km/sec 

Satellite  speed 

8  km/sec 

Electron  gyro  radius 

2.5  cm 

Ion  gyro  radius 

430  cm 

Electron  (Ion)  density 

10 11/m3 

Ion  gyro  radius 

0.74  cm 

Table  1:  Physical  parameters 


tion  at  a  boundary  point.  In  order  to  treat  the 
boundary,  we  assume  Maxwellian  velocity  distribu¬ 
tion1  for  electrons  at  the  boundary  and  translating 
Maxwellian  velocity  function  for  ions  at  the  far  up¬ 
stream.  First,  based  on  the  assumption  that  ions 
are  Maxwellian  in  the  far  upstream  region  and  are 
accelerated/decelerated  one  dimensionally,  we  show 
the  ion  distribution  function  at  the  boundary  as  a 
function  of  the  local  potential.  From  this  function, 
we  calculate  the  density,  which  is  required  for  the 
quasi-neutral  condition,  and  the  flux,  which  is  re¬ 
quired  to  calculate  the  number  of  ions  replenished 
into  the  computational  domain  per  timestep.  Elec¬ 
trons  are  assumed  to  have  a  Maxwellian  distribution 
at  the  boundary  due  to  the  fact  that  the  electron 
thermal  velocity  is  much  faster  than  the  satellite  or¬ 
bital  speed. 

The  parameters  used  in  this  computation  are  typ¬ 
ically  as  shown  in  Table  1. 


Ions 

Due  to  their  large  mass,  ions  are  very  slow  com¬ 
pared  to  the  satellite  speed.  This  fact  leads  us  to 
the  assumption  that  ions  are  Maxwellian  only  in  the 
far  upstream  region  and  that  each  ion  is  decelerated 
one-dimensionally  by  the  potential  field  created  by 
the  tether.  The  latter  assumption  may  be  verified 
by  the  fact  that  in  the  frame  moving  with  a  tether, 
the  major  velocity  component  of  an  ion  is  in  the  di¬ 
rection  of  the  tether  velocity,  but  of  opposite  sign. 
Even  if  the  ion  is  accelerated  in  a  random  direction, 
the  major  component  of  the  ion  velocity  is  still  the 
same.  Therefore  the  ion  trajectory  remains  almost 
straight,  leading  to  the  assumption. 

The  ion  density  is  obtained  as  follows.  Assum¬ 
ing  a  Maxwellian  distribution  for  ions  in  the  plasma 
frame,  we  have  an  ion  distribution  function  in  the 


1  Flowing  effects  on  electrons  are  discussed  and  formulated 
in  Appendix. 


far  upstream  region  as 


(u>x  ~  U )2 


+  WJ  +  wj 


} 


where  is  ion  density  at  infinity,  m*  ion  mass, 
Ti  ion  temperature,  wx,wy  and  wz,  x-,  y-  and  z - 
components  of  a  particle  velocity  respectively,  U 
the  tether  velocity  (in  the  x-directionj  and  vn  = 

KTi/rrii.  From  the  energy  conservation  between 
infinity  (the  far  upstream  region)  and  a  boundary 
point,  we  have 

^rmwl  =  ul  +  Ze<j>h  (2) 


Wy  ~  Wy 

WZ  ~  Wz 


(3) 

(4) 


where  wxis  an  x-component  of  a  particle  velocity  at 

a  boundary  point  and  is  a  boundary  potential. 
Substituting  (2)  into  (1),  we  get  the  ion  distribution 
function  at  the  boundary  as 


fion,b  —  £»  exp 


(5) 


where  &  =  n (m1/,27r/v711)3'/2.  Hereafter  we  omit 
the  hat  to  denote  the  quantities  at  the  boundary. 
From  this  function,  we  obtain  density  and  flux  of  in¬ 
coming  ions  as  a  function  of  the  boundary  potential. 
The  density  is  given  by 


n\n  =  f  Ai  exp 

Jo 


1 'Ti 


dwx 


(6) 


for  (fib  >  0,  and 

/oo 

A*  exp 

/  2  e<fih 

V" 


vTi 


dwx  (7) 


for  (fib  <  0,  where  the  integral  over  wy  and  wz  has 

been  carried  out,  and  Ai  =  noo  \fmi/2'KKTi.  The 
flux  is  then  given  as 


nn = f 

Jo 


XiWx  exp 


(v^g+^-tQ2! 


uTi 


dwx  (8) 


for  (fib  >  0  and 


pin 

**  i 


A iwx  exp 


dwx  (9) 


for  (fib  <  0.  In  consequence  of  the  approximation, 
we  assume  no  ions  coming  into  the  computational 
domain  from  downstream.  Therefore,  we  have 

njn  =  0  (10) 

T\n  =  0  (11) 


for  the  downstream  boundary. 


Computational  outside  boundary 


Figure  1:  Geometry  of  Computation 


Electrons 


From  Table  1  we  see  that  electrons  have  a  high  mo¬ 
bility  so  that  they  reach  the  thermal  equilibrium 
quickly.  Assuming  a  Maxwellian  distribution  for  in¬ 
coming  electrons  at  the  boundary,  we  have 


L,b  =  exp 


w*+w%  +  w* 


VTc 


exp 


(12) 


where  £e  =  ^oo  ( me/2nKTe )3^  .  When  06  <  0,  an 
electron  has  been  decelerated  by  the  time  it  reaches 
the  boundary.  Thus  the  minimum  speed  of  an  elec¬ 
tron  entering  into  the  computational  domain  is  still 

zero,  otherwise  it  is  \  To  facilitate  the  calcu- 

7  v 

lation  of  the  local  density  and  flux  of  incoming  elec¬ 
trons,  we  change  variables  from  ( wx ,  wy)  to  (wn ,  wt ) 
where  wn  is  the  normal  velocity  component  to  the 
boundary  surface,  wt  the  tangent  component  , 


When  06  >  0,  all  electrons  are  accelerated  including 
a  particle  corresponding  to  w  =  0  at  infinity  in  the 
plasma  frame.  This  particle  is  also  accelerated  by 
the  positive  potential,  raising  the  minimum  velocity 

of  incoming  electrons  up  to  Therefore  elec¬ 

trons  slower  than  this  velocity  are  prohibited  from 
entering  the  computational  domain.  The  density  of 
incoming  electrons  thus  becomes 

///  “9> 

Vexcept  for  /p  ^  } 


where 

V  :  {— oo  <  Wn  <  0;  —  oo  <  wt  <  oo;  —oo  <  wz  <  oo}  (20) 


Electron  flux  is  calculated  likewise.  For  06  <  0, 
we  have 


/oo  roc  r0 

I  I  -VJnfejdlVn  dwt  dwz 
-oo  J  — oo  j  —  oo 

■($) 


nooVTe 

2y/ir 


•exp 


(23) 

(24) 


And  for  06  >  0,  we  have 


<*> 

tr^r-T  V  1  e  * 


Vexcept  for  Jp 

After  some  algebra,  we  obtain 

f  =  =gr(t+1) 


(26) 


Quasineutrality  Condition 


wx  =  wn  cos  6  —  wt  sin  6  (13) 

Wy  =  wn  sin  0  +  wt  cos  e  (i4) 


where  0  is  the  angle  between  the  normal  of  the  sur¬ 
face  and  the  plasma  flow  (in  the  positive  x-direction) 
(See  Figure  1). 

Then  we  have  a  Maxwellian  distribution  function 
as 


/e,6 


we  exp 


rnc{w*  +w?  +  mf} 
2k  Te 


where 


(15) 


(16) 


For  06  <  0,  the  density  of  incoming  electrons  at 
the  boundary  is  given  by 


n 


fc^dwn  dwt  dwz 


(17) 

(18) 


As  shown  above,  densities  due  to  incoming  particles 
are  given  as  a  function  of  the  local  boundary  poten¬ 
tial  06.  In  [1],  densities  due  to  outgoing  particles  are 
calculated  numerically  by 


e,i  v--\fc  -»  t,  n 
2L,i=i  Windtb 

where  k  is  the  number  of  particles  going  out  of  the 
computational  domain  at  a  local  boundary  point  (in 
computaion,  it  is  a  cell),  W{  the  particle  velocity  go¬ 
ing  out,  n  the  normal  vector  of  the  boundary  surface, 
5  the  surface  area  ,  all  of  which  are  calculated  at  the 
same  boundary  point.  The  quasi-neutral  condition 
is  applied  by  equating  the  electric  charge  density  of 
outgoing  particles  and  incoming  particles  to  zero, 
Densities  of  incoming  particles  are  given  as  a  func¬ 
tion  of  a  local  boundary  potential  06. 

n°ut  —  n°ut  =  ni?(<pb)  —  n\n(4>b)  (28) 

The  boundary  potential,  06,  is  solved  numerically 
from  (28),  and  used  as  a  boundary  condition  in  a 
Poisson  solver. 


Results 

In  this  section,  we  show  the  results  from  our  simula¬ 
tions.  First,  we  consider  the  numerical  convergence 
and  stability.  Figure  2  shows  the  sum  of  total  en¬ 
ergies  of  all  particles  in  the  computational  domain. 
In  this  figure,  we  see  the  convergence  of  overall  to¬ 
tal  energy  of  particles  and  the  energy  increase  re¬ 
ported  in  [3]  is  not  identified.  This  is  because,  in 
our  simulations,  particles  refresh  their  energy  once 
they  either  are  absorbed  by  a  tether  or  go  out  of  the 
computational  domain,  and  then  are  recycled  for  the 
computation.  Therefore  particles  do  not  accumulate 
a  numerical  error  throughout  the  simulation. 

Figure  3  shows  the  number  of  particles  which  re¬ 
side  in  the  computational  domain.  It  also  indicates 
a  steady  state  and  does  not  show  an  increase  in  time 
as  reported  before  [[5]]. 

Figure  4  shows  typical  trajectories  of  electrons. 
The  effect  of  the  magnetic  field  can  be  seen  in  the 
ram  region  where  an  electron  is  trapped  in  the  field 
line.  In  and  neax  the  wake,  due  to  the  depletion  of 
ions  in  the  wake  (Figure  5),  the  electric  potential 
becomes  negative  with  respect  to  infinity  (Figure  6) 
so  that  most  of  the  electrons  are  reflected  by  this 
potential  bump.  In  the  immediate  vicinity  of  the 
tether,  electric  field  effects  become  dominant  over 
magnetic  field  effects.  Thus  the  magnetization  effect 
on  electrons  becomes  very  small. 

Due  to  the  high  positive  potential  on  the  tether, 
ions  are  decelerated  and  deflected  from  its  quasi  one¬ 
dimensional  trajectory,  giving  rise  to  a  narrow  layer 
in  the  shape  of  a  bow  shock  in  front  of  the  tether, 
where  ion  density  is  high  because  ion  trajectories  are 
tangent  to  it  (an  “acoustic”  surface)  (Figure  5).  The 
potential  at  the  ion  stagnation  point  is  equal  to  the 
ion  ram  energy  (5eV). 

In  order  to  maintain  the  quasi-neutrality,  elec¬ 
trons  are  attracted  to  the  ion  “bow  shock”  (Figure 
7).  In  Figure  8,  right  in  front  of  the  stagnation  point, 
quasi-neutral  electric  charge  density  is  recognized. 
To  attract  electrons  along  that  magnetic  field  line,  a 
positive  potential  “wing”  spreads  out  along  this  field 
line,  since  the  source  of  electrons  is  only  from  this 
direction,  whereas  electrons  residing  in  front  of  this 
potential  “wing”  are  trapped  by  the  magnetic  field, 
and  can  not  contribute  to  the  quasi-neutrality  near 
the  tether. 

Finally  Figure  9  shows  the  collected  current  and 
OML  current.  The  collected  current  turns  out  to 
be  twice  as  much  as  the  OML  current.  This  cur¬ 
rent  collection  may  be  partially  attributed  to  the 
fact  that,  although  the  geometry  is  2-dimensional  (a 
long  cylindrical  tether),  particles  are  accelerated  3- 
dimensionally  through  the  magnetic  gyration.  That 
is,  a  particle  which  has  an  in-plane  velocity  equal  to 
zero  at  infinity  can  have  an  in-plane  velocity  faster 


than  y  at  a  boundary,  whereas  without  the 
gyration,  this  is  the  velocity  it  must  have  at  the 
boundary.  This  leads  to  more  one-sided  flux  into 
the  domain,  hence  more  current  collection.  However 
further  analysis  should  be  performed  to  verify  the 
correctness  of  this  current  collection  result. 

In  conclusion,  we  present  in  this  paper  a  PIC  tech¬ 
nique  to  treat  the  boundary  for  flowing  magnetized 
plasmas.  This  method  shows  numerical  convergence 
and  stability.  And  as  a  preliminary  result,  we  obtain 
a  current  collection  which  is  twice  as  much  as  OML 
current. 
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Appendix:  Flow  effects  on  elec¬ 
trons 


In  this  paper,  we  assumed  a  Maxwellian  distri¬ 
bution  for  electrons  at  the  boundary.  When 
U  —  0  this  assumption  becomes  exact,  since  the 
Maxwellian  distribution  function  is  in  the  form  of 

exp  -^§0  where  E  =  \ mew 2  —  ecp.  Therefore  un¬ 
der  Maxwellian  collisionless  conditions,  this  is  the 
exact  solution  at  the  boundary.  However  for  U  ^  0, 
the  distribution  function  given  by 


fc,b  =  Ceexp 


(wx  -  U)2  +w%  +  w2 


(29) 


, where  £e  =  rioo  (me/27r/cTe)3/^2,  is  in  the  form  of 
exp  Then  the  wx  term  produces  an 

error  as  the  same  distribution  function  (29)  is  as¬ 
sumed  at  the  boundary  as  at  infinity.  In  this  ap¬ 
pendix,  the  simulation  using  the  shifted  Maxwellian 
distribution  (29)  at  the  boundary  is  shown.  And 
the  result  is  compared  with  the  case  of  Maxwellian 
distribution  (12)  and  the  error  is  discussed. 

As  we  derived  the  density  of  incoming  electrons 
from  (12)  to  (18)  and  (22)  and  the  flux  from  (12)  to 
(24)  and  (26),  we  likewise  derive  the  density  and  the 
flux  from  (29)  with 

V  :  {— oo  <  Wn  <  0;  —  oo  <  wt  <  oo;  — oo  <  wz  <  oo}  (30) 


Ip  :  (wn  —  U  cos  0)2  +  (tut  +  U  sin  0)2  -f  w2  <  — — 


me 


(31) 


The  result  is  as  follows;  The  density  is  for  U  cos  9  > 
\j2e4>fme  and  for  <frb  <  0 

exp  (S)  erfc(£ cos*)  (32) 
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And  the  flux  is  for  U  cos  0  >  ^2e<f)/me  and  for 
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Figure  2:  Total  energy  of  particles  in  a  domain 
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and  for  | U  cos0|  <  y/2e(j)/me 
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Figure  3:  Number  of  particles  in  a  domain 
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With  these  equations,  the  same  simulation  is  per¬ 
formed.  The  result  shows  the  current  collection 
which  is  20%  more  them  the  case  of  U  =  0  (Figure 
10).  From  this  we  can  deduce  that  the  error  of  using 
the  non-shifted  Maxwellian  distribution  function  at 
the  boundary  is  of  the  same  order. 
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Abstract 

This  thesis  carries  out  both  technical  and  policy  analyses  on  a  progression  of 
topics.  These  topics  are  cost  modeling  of  distributed  satellite  systems, 
determination  of  propulsion  system  requirements  for  satellite  clusters  and 
swarms,  and  an  analysis  of  ion  micropropulsion  systems  for  use  in  swarms  of 
microsatellites.  The  total  cost  over  a  10  year  mission  life  is  calculated  for 
configurations  of  the  NPOESS  mission  in  which  the  primary  instruments  were 
distributed  among  three  smaller  satellites.  This  increased  reliability  of  the 
distributed  configuration  substantially  reduces  the  number  of  ground  and  on- 
orbit  spares  required.  As  a  result,  mission  costs  are  reduced  compared  to  a 
single  large  satellite  configuration.  The  relative  positions  of  satellites  in  a 
cluster  are  altered  by  “tidal”  accelerations  which  are  a  function  of  the  cluster 
baseline  and  orbit  altitude.  Near  continuous  thrusting  by  a  propulsive  system 
was  found  necessary  to  maintain  the  relative  positions  of  the  satellites  within 
allowable  tolerances.  Satellite  mass,  volume,  and  power  constraints  limit 
reasonable  cluster  baselines  to  approximately  30  m,  300  m,  and  5000  m  at 
1000  km,  10,000  km,  and  GEO  altitudes  respectively.  To  maintain  these  cluster 
baselines,  the  propulsive  system  must  operate  at  specific  impulses  and 
efficiencies  consistent  with  those  of  ion  engines.  The  performance  of  the  linear 
ion  microthruster  concept  is  examined  using  Brophy’s  model  to  predict  energy 
costs  per  beam  ion.  Efficiency  of  the  linear  ion  microthruster  was  calculated  to 
be  in  the  range  of  15%.  At  this  efficiency,  the  linear  ion  microthruster  is  not  able 
to  perform  station  keeping  requirements  for  cluster  missions. 
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Chapter  1 

Introduction 


1.1  BACKGROUND 

This  thesis  addresses  the  issue  of  in-orbit  propulsion  for  micro-satellites  and 
distributed  satellite  systems  from  both  a  technical  and  policy  perspective.  The 
following  questions  will  be  answered: 

•  Are  there  satellite  missions  in  which  it  makes  sense,  from  a  cost 
perspective,  to  distribute  functionality  across  several  smaller 
satellites?  What  preliminary  design  considerations  should  satellite 
manufacturers  analyze  to  benefit  from  distribution? 

•  If  a  group  of  small  satellites  were  to  orbit  in  a  local  cluster,  what  are 
the  propulsion  requirements  necessary  to  maintain  that  formation? 

•  Is  it  feasible  to  scale  down  a  traditional  ion  engine  so  that  it  may  be 
used  to  maintain  cluster  formation  or  as  the  main  propulsion  for  a 
swarm  of  micro-satellites? 

This  chapter  reviews  work  performed  by  others  on  these  topics  and  provides 
background  on  the  analysis  carried  out  in  this  thesis. 

1.1.1  Cost  Modeling  of  Distributed  Satellite  Systems 

There  are  currently  two  satellite  systems  operated  by  U.S.  agencies  to  provide 
atmospheric  data  for  weather  forecasting  -  the  DMSP  system,  operated  by  the 
Department  of  Defense,  and  the  POES  system,  operated  by  the  Department  of 
Commerce.  Both  systems  are  approaching  the  end  of  their  operational  life,  and 
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the  two  departments  have  issued  a  common  Request  for  Proposal  [1]  for  the 
development  of  the  next  generation  weather  satellite  system.  Its  name  is 
NPOESS  (National  Polar-Orbiting  Operational  Environmental  Satellite 
System). 


The  mission  of  the  NPOESS  system  is  to  measure  and  observe  atmospheric 
and  space  environment  data  to  provide  accurate,  reliable,  and  up-to-date 
Environmental  Data  Records  (EDRs)  to  central  ground  stations  and  field  users 
spread  around  the  world.  The  NPOESS  requirements  specify  that  the 
instruments  should  fly  in  three  sun-synchronous  orbital  planes  with  nodal 
crossing  times  of  5:30  am,  9:30  am  and  1:30  pm. 


The  measurement  of  the  EDRs  specified  in  the  RFP  require  twelve  instruments 
on  each  orbital  plane,  three  of  which  are  considered  critical.  A  satellite  carrying 
a  primary  instrument  must  be,  by  definition,  replaced  if  that  instrument  fails. 
Distributing  the  primary  instruments  across  several  smaller  satellites  may 
increase  the  reliability  of  the  satellites  and  thus  the  system.  This  increased 
reliability  may  reduce  the  number  of  spares  required  thereby  reducing  the  cost 
of  the  system.  The  three  primary  instruments  and  the  distribution  configurations 
examined  in  this  thesis  are  shown  in  Figure  1-1. 


1  satellite  per  plane  1  satellite  per  plane 

1  set  of  critical  instruments  2  set  of  critical  instruments 

per  plane  per  plane 

^  ©  3sats/lset  .. 

tvTTFsl  Irvndl  frr^ol 

3  satellites  per  plane 
1  set  of  critical  instruments  per  plane 


Figure  1-1  Cluster  Configurations  Analyzed 
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1.1.2  Analysis  of  Cluster  Missions 

Future  satellite  missions  may  be  composed  of  several  satellites  flying  in 
formation  rather  than  a  single  larger  satellite.  Clusters  of  satellites  may  be 
employed  to  increase  the  reliability  of  the  system,  form  a  large  sparse  aperture, 
or  simply  to  provide  greater  coverage  of  an  area.  In  each  case,  a  propulsion 
system  will  be  required  to  maintain  the  formation  in  addition  to  normal  station¬ 
keeping. 

Jansen  [2]  reviewed  the  accelerations  experienced  by  satellites  orbiting  in  a 
formation.  In  particular,  he  examined  phase  control  at  all  altitudes,  creation  of 
“designer”  orbits  for  electrically  propelled  vehicles  as  well  as  formation  flying  of 
local  satellite  clusters.  Figure  1-2  depicts  a  cluster  of  satellites  forming  a  sparse 
aperture  for  high  resolution  imaging. 


Figure  1-2.  Satellite  Cluster  Formation 


Future  spacecraft  may  employ  microthrusters  for  missions  requiring  precise, 
low-thrust  firings.  These  missions  may  include  precision  station-keeping  of 
separated  spacecraft  forming  a  large  sparse  aperture,  control  of  large  flexible 
structures  such  as  deployable  antennas  or  solar  arrays,  or  as  the  main 
propulsion  system  for  microsatellites.  These  missions  indicate  a  need  to  greatly 
reduce  thrust,  size  and  power  from  those  of  current  engines  and  associated 
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hardware.  Additionally,  new  manufacturability  and  materials  problems 
associated  with  the  very  small  scales  must  be  overcome. 

1.1.3  Modeling  Performance  of  Micro  Ion  Engines 

Before  detailing  the  models  used  to  predict  the  performance  of  ion  engines,  it 
would  be  helpful  to  review  the  operation  of  the  thruster.  As  shown  in  Figure  1-3, 
ion  engines  are  cylindrical  chambers  with  diameters  typically  ranging  from  5  - 
30  cm.  Ring  magnets  are  placed  so  as  to  create  an  appropriately  shaped 
magnetic  field  within  the  chamber.  A  neutral  gas,  such  as  xenon,  is  fed  into  the 
discharge  chamber.  Electrons  are  injected  from  a  cathode  and  travel  within  the 
chamber  until  they  collide  with  a  neutral  atom  to  create  an  ion  or  until  they 
reach  an  anode.  The  electrons  tend  to  spiral  around  the  field  lines  increasing 
the  likelihood  that  a  collision  does  occur.  At  the  chamber  exit,  a  grid  system  is 
used  to  set  up  a  potential  difference.  Any  ions  in  the  neighborhood  of  the  grids 
will  be  accelerated  across  the  potential  difference  forming  an  ion  beam  and 
thus  thrust.  The  ion  beam  is  neutralized  by  an  exactly  opposite  charge  of 
electrons  collected  by  the  anode. 


v 

NET 


Figure  1-3.  Operation  of  Ion  Engines 
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Brophy  and  Wilbur  [3]  developed  a  model  to  predict  the  performance  of  ring 
cusp  ion  thrusters.  The  model  is  based  on  conservation  of  mass,  charge,  and 
energy.  The  result  is  a  single  algebraic  equation  formulated  in  terms  of  the 
average  energy  expended  in  producing  ions  in  the  discharge  plasma  and  the 
fraction  of  these  ions  extracted  in  the  beam.  The  average  plasma  ion  energy 
can  be  found  as  a  function  of  propellant  utilization  given  geometric  design 
parameters  of  the  thruster  and  the  type  of  propellant. 

Although  Brophy’s  model  is  very  simple,  determining  the  energy  expended  in 
producing  ions  and  the  fraction  of  these  ions  that  reach  the  beam  is  not  as 
simple.  Arakawa  [4]  developed  a  computer  algorithm  to  calculate  those 
parameters.  The  algorithm  calculates  the  magnetic  field  in  the  discharge 
chamber,  follows  the  path  of  electrons  along  the  field  lines,  determines  the 
distribution  of  ions,  and  then  inputs  this  information  into  Brophy’s  algebraic 
equation  to  determine  thruster  performance. 

The  ability  to  predict  ion  engine  performance  is  important  for  two  reasons.  First, 
it  allows  one  to  examine  how  performance  changes  as  the  size  of  the  chamber 
is  scaled  down  several  orders  of  magnitude.  It  also  allows  one  to  predict  the 
performance  of  new  concepts.  One  concept  is  for  a  micro-machined  ion 
propulsion  system  put  forth  by  members  of  the  Jet  Propulsion  Laboratory  [5]. 
The  thruster,  shown  in  Figure  1-4,  consists  of  several  linear  discharge 
chambers  situated  parallel  to  each  other.  Each  discharge  chamber  has 
dimensions  of  approximately  lOOjxm  X  300pm  x  10  cm.  The  walls  separating 
the  discharge  chambers  are  built  up  from  alternating  layers  of  conducting  and 
insulating  materials.  Contained  within  these  walls  is  the  ion  accelerator  system 
which  includes  the  screen  grid,  and  the  ion  accelerator  grid. 
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Figure  1-4.  Linear  Ion  Microthruster  Concept 

1.2  PREVIEW  OF  CHAPTERS 

•  Chapter  2  -  Cost  and  reliability  models  are  developed  to  analyze  the 
effects  of  distribution  on  the  National  Polar  Orbiting  Environmental 
Satellite  System  (NPOESS).  Life  cycle  costs  are  calculated  for  three 
instrument  distribution  configurations.  Critical  technologies  necessary 
for  deploying  distributed  satellite  systems  are  identified. 
Considerations  satellite  manufacturers  would  need  to  follow  in 
preliminary  design  trades  between  traditional  large  satellites  and 
distributed  satellites  are  presented. 

•  Chapter  3  -  An  analysis  of  the  station-keeping  requirements  of  a 
distributed  satellite  formation  known  as  a  local  cluster  is  carried  out. 
First,  a  determination  of  the  tolerances  to  which  the  individual 
spacecraft  must  maintain  relative  position  is  made.  Next,  these 
tolerances  are  related  to  the  required  thrust  level,  specific  impulse, 
efficiency,  etc.,  of  the  station-keeping  thruster.  Having  examined  the 
characteristics  of  thrusters  required  by  distributed  satellite  systems,  the 
rest  of  this  thesis  examines  the  performance  of  ion  microthrusters. 

•  Chapter  4  -  The  derivation  of  the  model  used  to  predict  the  performance 
of  ion  thrusters  is  laid  out.  Brophy’s  model  is  derived  from  conservation 
of  mass,  energy,  and  charge  principles.  Fundamentals  of  Arakawa’s 
algorithm,  the  computer  code  which  implements  Brophy’s  model,  are 
described. 


16 


•  Chapter  5  -  Arakawa’s  algorithm  is  used  to  examine  the  effect  on 
performance  of  scaling  a  traditional  7cm  cylindrical  ion  engine  down 
three  orders  of  magnitude  Simulations  with  no  magnetic  field,  a 
constant  magnetic  field,  and  a  magnetic  field  varying  inversely 
proportional  with  chamber  size  are  run.  Factors  influencing 
performance  as  scale  is  reduced  are  identified. 

•  Chapter  6  -  Arakawa’s  algorithm  is  modified  to  analyze  a  linear  ion 
microthruster  concept.  The  linear  concept’s  design  is  optimized  with 
respect  to  the  number  of  chambers  and  chamber  geometry. 

•  Chapter  7  -  The  questions  addressed  and  answers  provided  by  this 
thesis  are  summarized.  Conclusions  from  both  a  policy  and  technical 
perspective  are  identified. 


1.3  SUMMARY 

The  purpose  of  this  chapter  has  been  two-  fold.  1)  To  provide  an  introduction  to 
the  questions  addressed  in  this  thesis  and  2)  to  outline  the  procedure  for  how 
these  questions  will  be  answered.  Background  information  on  the  NPOESS 
mission,  satellite  cluster  propulsion  requirements,  and  models  used  to  predict 
ion  engine  performance  was  discussed.  The  topics  of  the  chapters  in  this  thesis 
were  also  previewed. 
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Chapter  2 

To  Distribute  or  Not  To  Distribute: 
A  Policy  Analysis 


The  previous  chapter  provided  background  on  the  topics  to  be  discussed  in  this 
thesis.  This  chapter  presents  a  policy  analysis  of  the  NPOESS  mission 
deployed  as  a  distributed  satellite  system.  A  model  is  developed  to  compare 
the  costs  of  deploying  a  satellite  mission  as  a  distributed  system  versus  a 
traditional  single  satellite.  It  also  details  the  considerations  satellite 
manufacturers  must  take  into  account  when  deciding  whether  to  deploy 
satellites  in  a  distributed  system.  For  the  total  design  of  the  NPOESS  mission 
as  a  distributed  system,  see  the  MIT  design  report  [1].  A  greatly  extended 
review  of  tradeoffs  and  considerations  regarding  distributed  satellite  systems 
can  be  found  in  the  book  chapter  by  Shaw  and  Yashko  [2]. 

2.1  INTRODUCTION  TO  DISTRIBUTED  SYSTEMS 

There  are  many  reasons  why  a  distributed  architecture  is  well  suited  to  some 
space  applications.  Unfortunately  the  arguments  for  or  against  distribution  are 
fraught  with  subjectivity  and  firmly  entrenched  opinions.  It  currently  seems  that 
most  of  the  satellite  design  houses  in  the  country  are  internally  split  between  the 
proponents  and  opponents  of  distribution.  Each  camp  supports  one  side  of  the 
debate  vehemently  and  can  find  a  seemingly  endless  stream  of  supporting 
arguments  to  back  their  claims.  The  “radicals”  claim  that  the  development  of 
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large  constellations  of  small  satellites  leads  to  economies  of  scale  in 
manufacture  and  launch,  reducing  the  initial  operating  costs.  They  also 
expound  that  the  system  becomes  inherently  more  survivable  due  to  the  in-built 
redundancy.  Conversely,  the  “traditionalists”  debunk  these  arguments, 
reminding  everyone  that  you  can’t  escape  the  need  for  power  and  aperture  on 
orbit,  and  that  building  even  100  satellites  does  not  imply  significant  bulk¬ 
manufacturing  savings.  They  assert  that  the  lifetime  operating  costs  for  large 
constellations  will  far  outweigh  the  savings  incurred  during  construction  and 
launch. 

In  fact,  most  of  the  statements  made  by  both  sides  are  true,  but  only  when  taken 
in  context.  Clearly  a  distributed  architecture  is  not  the  panacea  for  all  space 
applications.  It  is  tempting  to  get  carried  away  with  the  wave  of  support  that  the 
proponents  of  distributed  systems  currently  enjoy.  Care  must  be  taken  to  curb 
this  blind  faith.  Also  best  avoided  is  the  naive,  but  commonplace  application  of 
largely  irrelevant  metaphors  supporting  the  adoption  of  distributed  systems;  the 
unerring  truth  that  ants  achieve  remarkable  success  as  a  collective  is  really  not 
an  issue  in  satellite  system  engineering! 


2.2  CONSIDERATIONS  FOR  IMPLEMENTING  DISTRIBUTED 
SATELLITE  SYSTEMS 

There  are  some  factors  that  are  critical  to  the  design  of  a  distributed  architecture 
that  were  irrelevant  to  the  design  of  traditional  systems.  Depending  on  the 
application,  these  issues  may  be  minor  hurdles,  or  could  be  so  prohibitive  that 
the  adoption  of  a  distributed  architecture  is  unsuitable  or  impossible.  Some  of 
the  important  considerations,  characteristic  of  all  distributed  architectures,  and 
particular  to  small-  and  microsatellite  designs  are  presented  here. 

The  distribution  of  system  functionality  among  separate  satellites  means  that  the 
system  is  essentially  transformed  into  a  modular  information  processing 
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network.  One  beneficial  aspect  of  modularization  comes  from  an  improved 
fault-tolerance.  System  reliability  is  by  nature  hierarchical  in  that  the  correct 
functioning  of  the  total  system  depends  on  the  availability  of  each  of  the 
subsystems  and  modules  of  which  the  system  is  composed.  Early  reliability 
studies  [3]  showed  that  the  overall  system  reliability  was  increased  both  by 
applying  protective  redundancy  at  as  a  low  a  level  in  the  system  hierarchy  as 
was  economically  and  technically  feasible,  and  by  the  functional  separation  of 
subsystems  into  modules  with  well-defined  interfaces  at  which  malfunctions  can 
be  readily  detected  and  contained.  Clearly,  subdividing  the  system  into  low- 
level  redundant  modules  leads  to  a  multiplication  of  hardware  resources  and 
associated  costs.  However,  the  impact  of  improved  reliability  over  the  lifetime  of 
the  system  can  outweigh  these  extra  initial  costs. 

There  are  additional  factors  supporting  modularization  that  are  specific  to 
satellite  systems.  The  baseline  costs  associated  with  a  system  of  small 
satellites  may  be  smaller  than  for  a  larger  satellite  design.  Of  even  greater 
impact  is  the  lower  replacement  costs  required  to  compensate  for  failure.  A 
modular  system  benefits  not  only  because  a  smaller  replacement  component 
has  to  be  constructed,  but  also  because  of  the  huge  savings  in  its  deployment. 

All  of  these  factors  suggest  that  a  system  should  be  separated  into  modules  that 
are  as  small  as  possible.  However,  there  are  some  distinct  disadvantages  of 
low-level  modularization  that  must  be  considered.  The  most  important  of  these 
are  the  costs  and  low  reliability  associated  with  complexity. 


2.2.1  Complexity 

The  complexity  of  a  system  is  well-understood  to  drive  the  development  costs 
and  can  significantly  impact  system  reliability.  In  many  cases,  complexity  leads 
to  poor  reliability  as  a  direct  result  of  the  increased  difficulty  of  system  analyses; 
failure  modes  were  missed  or  unappreciated  during  the  design  process.  For  a 
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system  with  a  high  degree  of  modularity,  these  problems  can  offset  all  of 
benefits  discussed  above. 

Although  each  satellite  in  a  distributed  system  might  be  less  complex,  being 
smaller  and  having  lower  functionality,  the  overall  complexity  of  the  system  is 
greatly  increased.  The  actual  level  of  complexity  exhibited  by  a  system  is 
difficult  to  quantify.  Generally,  however,  it  is  accepted  that  the  complexity  is 
directly  related  to  the  number  of  interfaces  between  the  components  of  the 
system.  Although  the  actual  number  of  interfaces  in  any  system  is  architecture 
specific,  it  is  certainly  true  that  a  distributed  system  of  many  satellites  has  more 
interfaces  than  a  single  satellite  design.  Network  connectivity  constraints  mean 
that  the  number  of  interfaces  can  increase  geometrically  with  the  number  of 
satellites  in  a  distributed  architecture.  This  is  an  upper  bound;  systems 
featuring  satellites  operating  in  parallel  with  no  inter-satellite  communication 
(defined  as  collaborative  systems)  exhibit  linear  increases  in  interfaces  with 
satellites.  The  complexity  of  a  distributed  system  is  therefore  very  sensitive  to 
the  number  and  connectivity  of  the  separate  modules. 

2.3  MODELING  DISTRIBUTED  SATELLITE  SYSTEMS 

In  1996,  a  request  for  proposals  (RFP)  was  issued  for  the  next  generation 
National  Polar  Orbiting  Environmental  Satellite  System  (NPOESS)  [4].  This 
provided  an  opportunity  to  examine  the  merits  of  deploying  NPOESS  with 
functionality  distributed  across  several  smaller  satellites.  The  nodal  crossing 
time  requirements  of  the  RFP  dictate  that  the  NPOESS  space  segment  consist 
of  three  orbital  planes.  There  is  no  a  priori  reason,  other  than  co-registration 
requirements,  that  all  sensors  in  the  plane  must  be  located  on  a  single  satellite. 
In  fact,  locating  all  sensors  on  a  single  satellite  reduces  the  reliability  of  the 
satellite  since  failure  of  any  of  the  several  primary  sensors  would  require 
replacement  of  the  entire  satellite.  It  may  be  possible  then,  to  increase  the 
reliability  and  reduce  the  lifetime  cost  of  the  NPOESS  space  segment  by 
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distributing  the  sensors  across  several  smaller  satellites  orbiting  as  a  cluster. 
This  chapter  describes  the  model  developed  to  examine  several  candidate 
cluster  configurations  to  determine  the  effect  that  distributing  the  primary 
sensors  has  on  total  lifetime  costs. 

2.3.1  Model  Overview 

The  methodology  of  the  reliability  /  cost  model  is  shown  in  Figure  2-1 .  The 
number  of  satellites  in  each  plane  along  with  the  types  of  sensors  on  each 
satellite  are  chosen  as  primary  inputs  to  the  model.  In  addition,  the  bus  design 
life,  bus  reliability,  sensor  reliability,  mission  life,  learning  curve  percent,  and  a 
term  known  as  the  level  of  similarity  (LOS)  must  be  included.  For  configurations 
where  instruments  are  distributed  across  several  spacecraft,  the  buses  would 
be  similar  but  not  identical.  The  level  of  similarity,  ranging  from  0  to  1,  is  a 
measure  of  how  similar  the  buses  are  so  that  costs  may  be  shared  where 
appropriate.  Probability  of  bus  failure  is  assumed  to  be  100%  at  the  end  of  the 
bus  design  life  due  to  exhaustion  of  expendables,  failure  of  batteries,  or 
insufficient  power  due  to  solar  array  degradation. 


r#  Number  ot  Satellites  > 
in  formation 
•  Types  of  sensors  on 
each  satellite  type 
(ie.  mass,  power,  reliability) 


lifetime  Costs  faetuimg  In: - \ 

-  Learning  Curve  from 
multiple  satellite  production. 

-  Discounting  due  to  satellite 
replacements  in  the  future. 


r  Reliability  of  each 
type  of  satellite  in 
cluster. 


M;Ks,  power,  size  of 
each  type  of  satellite 
in  cluster. 


fRDl&band  1HJ 
cost  of  each 
satellite  type. 


Failure  densities 
of  each  satellite  type. 


Total  number  of  satellites 
of  each  type  required  over 
10  yr.  mission. 

Satellite  replacement 
schedule. 


Figure  2-1  Cost  Model  Methodology 


23 


Bus  subsystem  characteristics  are  estimated  from  payload  mass  and  power 
using  estimates  based  on  a  regression  of  data  from  15  previous  Department  of 
Defense  communication,  navigation,  and  remote  sensing  missions.  With  the 
bus  and  payload  information,  the  reliability  of  each  spacecraft  type  can  be 
calculated.  Failure  rates,  and  therefore  an  expected  deployment  schedule,  for 
each  satellite  type  follow  from  the  spacecraft  reliability.  Parametric  costs 
models  can  then  be  used  to  estimate  the  cost  of  the  sensors  and  the  spacecraft 
bus  based  upon  the  estimated  subsystem  mass  and  power  requirements. 
Knowing  the  number  of  satellites  deployed  and  their  expected  launch  date,  the 
total  cost  over  the  system’s  mission  life  can  be  calculated  factoring  in  learning 
curve  effects  and  discounting  to  current  year  dollars.  The  details  of  this  model 
are  described  in  subsequent  sections  of  this  chapter. 

2.3.2  Payload  Cost  Model 

An  estimate  of  the  cost  of  each  sensor  was  calculated  using  the  Goddard  Space 
Flight  Center  (GSFC)  Multi-Variable  Instrument  Cost  Model.  The  GSFC 
instrument  cost  model  generates  a  cost  estimate  to  develop  and  produce  the 
first  unit  from  a  regression  of  sensor  type,  mass,  power  and  data  rate  for  186 
existing  instruments.  An  additional  factor  is  included  to  account  for  the  level  of 
technology  in  the  year  the  instrument  was  developed.  The  cost  to  procure  the 
sensor  is  assumed  to  be  50%  of  the  prototype  cost.  The  parametric  relationship 
used  to  estimate  the  prototype  cost  of  the  instrument  (FY97$M)  is 

CO5r  =  0.576((0.453MASS)a421)(pWRa,77)((i7?-1960)~0113)(Di?r0048)(FAM102®)(CLSa928)  (2-1) 

where  MASS=  instrument  mass  (kg) 

PWR  =  instrument  power  (W) 

DRT  =  instrument  data  rate  (kbps) 

YR  =  year  launched 
FAM  =  mission  family 
CLS  =  instrument  class 

Valid  ranges  for  instrument  mass,  power,  data  rate  and  launch  year  are  2.2  to 
6847.5  kg,  0.6  to  2000  W,  0.008  to  85000  kbps,  and  1965  to  1989  respectively. 
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In  general,  it  is  not  advisable  to  extrapolate  more  than  25%  beyond  a 
parameter's  range  [5].  It  should  be  mentioned  that  the  NPOESS  sensors,  which 
were  assumed  to  consist  of  1997  technology,  represent  an  extrapolation  of  29% 
beyond  the  range.  Values  for  the  instrument  class,  CLS,  and  mission  family, 
FAM,  are  listed  in  Table  2-1.  All  sensors  used  onboard  NPOESS  satellites 
were  assigned  a  Mission  Family  designation  of  “Normal  (Earth  resource  sat)”. 
Table  2-2  lists  the  Instrument  Class  and  designation  of  each  of  the  sensors.  For 
suites  consisting  of  sensors  belonging  to  more  than  one  instrument  class,  an 
average  CLS  value  was  used.  GPSOS  and  S&R,  which  do  not  specifically 
belong  to  one  of  the  GFSC  model  instrument  classes  were  assigned  to  an  F IF 
class  with  a  CLS  value  of  1 . 

Table  2-1  GSFC  Cost  Model  CLS  and  FAM  values 


Instrument  Class 

CLS 

Plasma  Probe 

1.030 

Magnetometer 

1.065 

Passive  Microwave 

1.182 

1.221 

Mass  Measurement 

1.270 

Electric  Field 

1.290 

Active  Microwave 

1.291 

Interferometer 

1.309 

EM— 

1.317 

Photometer 

1.490 

Laser 

1.630 

Radiometer 

Telescope 

2.300 

i;iMl=TAiAM.!LM 

4.890 

Mission  Family 

FAM 

Shuttle 

0.440 

Explorer  (Small  Free  Flier) 

1.125 

Normal  (Earth  resource  sat) 

1.720 

Planetary 

2.282 

. . . 

2.380 

Table  2-2  GSFC  Cost  Model  Instrument  Class  and  Designation 


Sensor  Name 

Instrument  Class 

VJIRS 

i maannHH 

primary 

CRiS 

Radiometer 

primary 

CMIS 

Passive  Microwave,  Active  Microwave 

GPSOS 

RF 

secondary 

S&R 

RF 

secondary 

SOBEDS 

Magnetometer,  Plasma  Probe,  Charge  detector 

secondary 

NCERES 

Radiometer 

secondary 

TOMS 

Spectrometer 

secondary 

Altimeter 

Radiometer 

secondary 

DIDM 

Electric  Field  Measurement,  Mass  Measurement, 
Charge  Detector 

secondary 
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NACRIM 

Radiometer 

secondary 

ABIS 

Radiometer 

secondary 

SUVPHO 

Photometer 

secondary 

ARGOS 

Radiometer 

secondary 

2.3.3  Satellite  Bus  Cost  Model 

Space  segment  costs  were  estimated  using  the  U.S.  Air  Force  (USAF) 
Unmanned  Spacecraft  Cost  Model  v5  [6].  This  model  utilizes  parametric  cost 
estimating  relationships  (CER)  derived  from  past  satellite  programs  based  on 
major  subsystem  characteristics  such  as  mass,  power,  and  data  rate.  These 
subsystem  parameters  were  estimated  using  average  values,  listed  in  Table  2- 
3,  of  15  previous  Department  of  Defense  communications,  navigation,  and 
remote  sensing  satellite  missions.  Payload  mass  and  power  are  known  given 
the  sensors  required  to  measure  the  environmental  data  records.  The  dry  mass 
of  the  spacecraft  bus  along  with  the  individual  subsystem  masses  can  then  be 
estimated. 


Table  2-3  Estimates  of  subsystem  masses 


Subsystem 

%  of  spacecraft 
dry  mass 

Payload 

28.0 

Structures 

21.0 

Thermal 

4.5 

Power 

30.0 

Telem,  Tracking,  Control 

4.5 

Attitude  Control 

6.0 

Propulsion  j 

6.0 

In  addition,  the  average  spacecraft  power  and  required  array  power  can  be 
estimated  from  the  total  payload  power  (Watts)  as  follows. 


P  = 

avg 


_  r  payload 


0.3 


(2-2) 


/»  =  1„33jP 

array  avg 


(2-3) 
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The  subsystems’  mass  and  power  can  be  used  to  estimate  the  subsystems’ 
research,  development,  test  and  evaluation  (RDT&E)  and  theoretical  first  unit 
(TFU)  cost  according  to  the  CER’s  shown  in  Table  2-4. 


Table  2-4  USAF  Unmanned  Spacecraft  RDT&E  and  TFU  CERs 


Cost  Component 

| 

ESe9m 

"HUTEE  CEE — 

(FY97$k) 

THT'CER - 

(FY97$k) 

BflfrVmKh— 1 — 

|  7  -  428 

mmamzm 

|&!QI!2L^M1 

TT&C 

|  Mass  (kg) 

Attitude  & 

Reaction  Control 

Rower  * 

M5BIBEE«asMM 

kLOC 

N/A  ] 

ESSE 

Aerospace 

Ground 

Equipment 

HUI&b  +  TFU  , 

Hardware  Costs 
($k) 

'23531  - 
285576 

n7a 

Program  Level 

Satellite  Hardware 
Cost  ($k) 

34012- 

267347 

■Ml 

n^im 

*  from  USAF  Unmanned  Spacecraft  Cost  Model  v.6 


The  space  segment  was  assumed  to  have  a  minimum  impact  on  existing 
ground  stations  regardless  of  the  number  of  satellites  in  the  configuration. 
Communications  subsystems  and  satellite  configuration  were  designed  to 
accommodate  existing  ground  stations  and  field  user  capabilities.  Major 
ground  stations  should  need  to  utilize  two  or  three  of  their  existing  antennas  to 
communicate  with  the  space  segment.  Moreover,  the  duration  of  satellite  data 
transmission  should  not  excessively  burden  ground  stations.  The  deployed  field 
users  require  only  a  single  standard  issue  antenna  for  any  of  the  proposed 
configurations.  Based  upon  these  assumptions  and  because  costs  vary  widely, 
the  cost  of  ground  stations  was  not  considered  in  the  configuration  trade. 

The  cost  to  launch  each  satellite  was  determined  by  examining  the  vehicles 
capable  of  being  launched  from  Vandenberg  Air  Force  Base  (VAFB)  listed  in 
Table  2-5.  The  wet  masses  of  the  satellites  were  compared  to  the  launch 
capacities  of  the  available  vehicles.  For  the  scenario  in  which  each  plane 
contains  a  cluster  of  more  than  one  satellite  type,  as  many  satellites  as  possible 
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were  placed  on  a  launch  vehicle  during  initial  deployment  of  the  system.  It  must 
be  kept  in  mind  that  the  final  launch  into  a  plane  may  not  need  to  carry  as  many 
satellites  as  the  previous  launches.  For  example,  if  a  plane  contains  a  cluster  of 
three  small  satellites  and  the  launch  vehicle  is  capable  of  carrying  two  satellites 
at  once,  then  the  second  launch  of  the  initial  deployment  into  that  plane  would 
carry  only  one  satellite.  Thus,  the  final  initial  deployment  launch  into  a  plane 
and  any  launch  to  replace  a  failed  satellite  may  occur  on  a  different  launch 
vehicle.  The  chosen  launch  vehicle  is  the  one  that  minimizes  launch  costs 
given  the  wet  masses  of  the  satellites. 


Table  2-5  Launch  Vehicle  Costs  and  Capacities 


Launch 

Vehicle 

Capacity  to 
SSPO  (kg) 

Cost  /  Launch  (97  M$) 

LMLV  II 

1210 

r  21 

Taurus  XL 

1150 

38 

Titan  II  S 

3028 

36 

Delta  II  7925 

3175 

60 

Titan  IV 

13364 

300 

2.3.4  Similarity  Factor 

The  RDT&E  cost  for  each  satellite  type  can  be  calculated  using  the  USAF 
unmanned  spacecraft  cost  model  described  in  Section  B-2.  In  the  case  where 
there  is  more  than  one  satellite  type,  a  separate  RDT&E  cost  would  be 
calculated  for  each  type.  Because  there  are  some  similarities  between  the 
buses,  each  satellite  type  would  not  have  to  bear  the  full  development  cost  of  a 
stand  alone  program.  The  level  of  similarity,  LOS,  ranging  from  0  to  1  indicates 
the  commonality  among  the  different  bus  types.  A  LOS  equal  to  1  represents 
identical  buses.  A  LOS  equal  to  0  indicates  that  no  RDT&E  costs  are  shared 
among  the  different  bus  types.  A  similarity  factor,  SF,  by  which  the  stand-alone 
RDT&E  cost  of  satellite  bus  type  b  is  multiplied  to  reflect  its  true  RDT&E  cost  in 
conjunction  with  other  similar  buses  under  development,  is  given  by: 


SFh=l- 


i__RD7^ 


(RDTEb  )(numbus)  J 


LOS 


(2-4) 
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where  numbus=  number  of  different  bus  types 

RDTEavg=  average  stand-alone  RDT&E  cost  of  all  bus  types 
RDTEb  =  stand-alone  RDT&E  cost  of  bus  type  b 

2.3.5  Learning  Curve  And  Discounting 

Two  additional  considerations  must  be  taken  into  account  when  determining  the 
costs  over  the  system’s  life.  The  learning  curve  is  a  method  to  account  for  the 
productivity  improvements  as  a  larger  number  of  satellites  are  produced. 
Included  in  this  concept  are  cost  reductions  due  to  economies  of  scale,  set-up 
time,  and  human  learning  as  the  number  of  satellites  produced  increases.  The 
cost  to  produce  the  next  unit  of  satellite  type  b  is  given  by 


pcos,b  =  TFUt[Nb  +(Nm  -Nt)LOS]°  -7TO,[(JV,  -l)+(AT_  -AT.jLOS]0  (2-5) 

where  Nb  =  number  of  satellite  bus  type  b  produced  including  current  unit 
Ntot  =  total  number  of  all  bus  types  produced  including  current  unit 
ln(100%/S) 

In  2 

S  =  slope  of  the  learning  curve  (%) 

TFU  =  theorectical  first  unit  cost 

The  LOS  term  is  included  to  allow  for  reductions  in  the  production  cost  of  a  bus 
type  due  to  the  production  of  other  similar,  yet  not  identical,  bus  types.  The 
learning  curve  slope,  S,  represents  the  percentage  reduction  in  cumulative 
average  production  cost  when  the  number  of  satellites  produced  is  doubled.  A 
95%  learning  curve  slope  was  applied  to  production  of  satellite  buses.  When 
calculating  the  expected  replacement  costs,  the  next  unit  is  represented  by  the 
probability  of  needing  to  build  a  satellite  bus  in  order  to  maintain  required 
number  of  spares. 


The  year  each  satellite  is  launched  can  be  determined  from  the  deployment 
schedule  discussed  in  the  next  section.  The  production  cost  of  each  satellite  is 
assumed  to  be  distributed  over  a  three  year  period  prior  to  the  launch  date.  20% 
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of  the  cost  is  incurred  two  years  before  launch,  30%  is  incurred  the  year  before, 
and  50%  is  incurred  during  the  year  of  the  launch  [7].  Furthermore,  the  cost  is 
discounted  to  constant  1997  dollars  by  multiplying  the  cost  in  each  year  by  the 
factor 


PV  = 


(1  +  d) 


n-1 


(2-6) 


where  n  =  the  year  the  cost  is  incurred  (relative  to  the  constant  dollar  year) 
d  =  discount  rate 


2.3.6  Replacement  Model 

The  failure  of  the  spacecraft,  and  therefore  its  replacement  schedule,  is  a 
function  of  the  reliability  of  the  spacecraft.  Because  a  spacecraft  must  be 
replaced  if  any  of  its  primary  instruments  fail,  this  model  assumes  that  the 
reliability  of  the  spacecraft  is  related  to  the  reliability  of  the  bus  and  primary 
sensors  as  follows: 


^s/c  ^ bus 


types 


i=l 


(2-7) 


where  =  reliability  of  spacecraft 
Rbus  =  reliability  of  bus 
R;  =  reliability  of  primary  instrument  i 
N;  =  number  of  primary  instruments  of  type  i  on  spacecraft 
types  =  number  of  different  types  of  primary  instruments  on  spacecraft 


It  is  required  that  a  99%  confidence  of  achieving  95%  availability  be  met  A 
Monte  Carlo  simulation  of  the  mission  was  used  to  calculate  failure  densities  at 
a  given  time.  These  failure  densities  are  then  used  to  determine  the  number  of 
spares  required  to  meet  the  confidence  and  availability  requirements.  For 
example,  if  the  probability  of  three  satellites  failing  in  the  time  it  would  take  to 
build  another  satellite  was  greater  than  1%,  then  three  spares  would  be 
required. 
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2.4  RESULTS  OF  CLUSTER  CONFIGURATION  ANALYSIS 

The  reliability/cost  model  is  used  to  determine  satellite  configurations  which 
best  meet  the  measures  of  effectiveness.  Figure  2-2  shows  the  three  candidate 
cluster  configurations  that  were  analyzed.  Three  primary  instruments,  VIIRS, 
CMIS,  and  CrIS,  are  required  for  the  mission.  Details  of  these  instruments  are 
discussed  in  reference  1.  Two  cluster  configurations  examined  consisted  of  the 
three  instruments  placed  either  all  on  a  single  satellite  (configuration  1)  or  each 
on  one  of  three  smaller  satellites  (configuration  3).  Another  configuration,  in 
which  an  additional  primary  sensor  is  added  to  the  single  satellite  for 
redundancy  (configuration  2),  was  also  examined. 


O  1  sat  /I  set 


@  1  sat  /2  sets 


fvnRsO  [cmisI]  fetus  [ 


VnRsIl  fcMisfl  fetus [ 


VnRsIl  fcMisfl  fetus  0 


1  satellite  per  plane 
1  set  of  critical  instruments 
per  plane 


©  3 


1  satellite  per  plane 

2  set  of  critical  instruments 

per  plane 


sats  /I  set 


Umi  fcMtsfl  fetus  ll 


3  satellites  per  plane 
1  set  of  critical  instruments  per  plane 


Figure  2-2  Cluster  Configurations  Analyzed 


For  configurations  consisting  of  more  than  one  satellite,  the  secondary 
instruments  were  distributed  among  the  satellites  in  such  a  way  as  to  balance 
payload  mass  and  power  requirements.  This  implies  a  fairly  standard  design 
and  a  high  level  of  similarity  amongst  the  buses.  Table  2-6  lists  the  sensors 
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along  with  their  respective  mass,  power  and  cost  characteristics  for  each 
configuration  shown  in  Figure  2-2. 


2.4.1  Cost  Breakdown  for  Cluster  Configurations 

Using  these  payload  values,  the  spacecraft  mass  and  power  characteristics  can 
be  found  using  the  relationships  described  in  Table  2-3  and  Equations 
(2-2)  through  (2-3).  These  spacecraft  mass  and  power  values  are  listed  in 
Table  2-7  for  each  configuration. 


Table  2-8  and  Table  2-9  list  the  subsystem  costs  for  the  1  sat  / 1  set,  1  sat  /  2  set 
and  3  sats  / 1  set  configurations.  These  values  were  estimated  using  the  CER’s 
in  Table  2-4.  The  total  RDT&E  cost  for  the  3  sats  / 1  set  configuration  must  be 
multiplied  by  the  similarity  factor  as  calculated  by  equation  (2-4). 


Table  2-6  Distribution  of  sensors  for  cluster  configurations 


Sensor 

iUBsaa 

1  sat/ 

1  sat/ 

3  sats  / 1  set 

Name 

mmm 

BS9 

Id!)] 

1  set 

2  set 

sat  1  sat  2  sat  3 

VII RS 

132 

215 

22.4 

V 

✓✓ 

✓ 

CMIS 

178 

208 

nm 

S 

V' 

CrIS 

68 

82 

7.1 

S 

ss 

✓ 

GPSOS 

9 

13 

1.2 

V 

s 

S 

S&R 

82 

88.8 

3.5 

S 

s 

✓ 

SOBEDS 

5 

5.6 

1.1 

S 

s 

✓ 

S 

✓ 

NCERES 

40 

35 

3.8 

s 

s 

S 

TOMS 

33 

25 

2.9 

V 

s 

V 

Altimeter 

53 

113 

5.3 

S 

s 

S 

NACRIM 

39 

35 

WEM\ 

V 

s 

s 

ABIS 

15 

7 

1.8 

</ 

V 

s 

SUVPHO 

8 

10 

1.4 

</ 

s 

ARGOS 

71 

62.4 

5.0 

✓ 

I  ✓ 

V 

Payload  Mass  (kg) 

733 

1110 

268 

263 

230 

Payload  Power  (W) 

900 

900 

358 

289 

290 

Payload  Cost  (FY97$M) 

68 

106 

32 

16 

24 

32 


Table  2-7  Spacecraft  mass  and  power  characteristics 


Characteristic 

1  sat/ 

1  set 

1  sat/ 

2  sets 

3  sats  / 1  set 
sat  1  sat  2  sat  3 

Dry  mass  (kg) 

2617 

3967 

955 

940 

821 

Wet  mass  (kg) 

2647 

3997 

985 

970 

851 

Avg.  Power  (W) 

2999 

2999 

1191 

963 

969 

Peak  Power  (W) 

3989 

3989 

1584 

1281 

1288 

BOL  Power  (W) 

5186 

5186 

2060 

1666 

1675 

Table  2-8  1  sat  / 1  set  and  1  sat  /  2  sets  subsystem  costs 


Configuration 
Cost  Component 

Payload 

Structured  hermal 

TT&C - 

Attitude  Determ. 

Reaction  Control 

Power 

Software 
Aeros  C3md  Equip 

Program  Level 
Launch  Ups  &  Orbital 
Support _ 


1  sat  / 1  set  H 

RDT&E  TFU - 1 

(FY97$M)  (FY97$M) 


"R77T 

EES 

EES 

TOT 

EES 

EES 

IT 

3BT 

E9T 

W 


~ss~ 

~ET 
1 OX 
SEE 
~SE~ 
■5T W 

~mr 

"R77T 

EOT 

TS~ 


1  sat  /  2  sets 


RDT&E 
(FY97$M) 
— WE~ 
"Tin — 
"  44.0 
48:5 
43.9 


— ITT 
(FY97 
“TOS 
— 91 

~~ m 

~1T‘ 

TTT 


33X 

XX 

03T 

93T 

W 


ITT 

TJ7A" 

Tf77T 

75T 

TIT 


Table  2-9  3  sats  / 1  set  subsystem  costs 
(RDT&E  costs  must  be  multiplied  by  similarity  factor) 


H; 

- Satellited - 1 

- Sate!! 

(FY97$M) 


(FY97$M) 


(FY97$M) 


(FY97$M) 


(FY97$M) 


fie; 


(F 


Table  2-10  lists  important  assumptions  inputted  to  the  reliability/cost  model. 


Table  2-10  Cost  Model  Assumptions 


Sensor  Reliability 

Mission  Life 

10  years  j 

Bus  Reliability 

HHHEEEHHH 

Launch  Vehicle  Reliability 

0.98 

Bus  Level  of  Similarity  (LOS) 

0.66 

Discount  Rate 

6% 

Year  of  Initial  Deployment 

2004 

2.4.2  Total  Mission  Costs 

The  total  costs  over  a  10  year  mission  life  were  calculated  for  each  of  the  three 
cluster  configurations.  As  shown  in  Figure  2-3,  the  costs  over  the  10  year 
period  are  broken  up  into  four  categories;  namely  RDT&E,  initial  deployment, 
required  spares,  and  expected  replacements.  RDT&E  cost  for  the  3  sats  /  1  set 
configuration  are  those  listed  in  Table  2-9  multiplied  by  the  similarity  factor. 
Initial  deployment  includes  the  development,  production,  and  launch  costs  for 
each  orbital  plane’s  original  complement  of  spacecraft.  The  number  of  required 
bus,  payload,  and  launch  vehicle  spares  is  derived  from  the  operations  model 
so  as  to  meet  required  mission  availability  (95%)  with  the  desired  confidence 
(99%).  Expected  replacements  also  flow  from  the  operations  model  and 
indicate  the  number  of  bus,  payload,  and  launch  vehicle  units  which  must  be 
produced  during  the  mission  to  maintain  the  required  number  of  spares. 

Figure  2-3  shows  that  the  initial  deployment  cost  is  least  expensive  for  the  1  sat 
/  1  set  configuration.  Adding  a  redundant  sensor  to  the  single  satellite 
configuration  greatly  increases  initial  deployment  cost  in  terms  of  larger  bus 
size,  additional  instruments,  and  more  expensive  launch  vehicles.  The  3  sat  / 1 
set  configuration,  although  being  launched  on  a  less  expensive  vehicle,  is 
slightly  more  expensive  than  the  1  sat  /  1  set  configuration  due  to  the 
duplication  of  bus  subsystems  (e.g.  power,  attitude  control,  propulsion  etc.)  and 
some  sensors  (e.g.  GPSOS)  in  each  orbital  plane. 
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Figure  2-3  Total  costs  over  10  year  mission  life 


Figure  2-3  also  shows  that  adding  a  redundant  sensor  increases  the  cost  as 
compared  to  configurations  with  a  single  primary  instrument.  The  slight 
decrease  in  the  failure  densities  as  a  result  of  redundancy  does  not  make  up  for 
the  expense  of  additional  sensors.  It  is  assumed  that  periods  of  one  year  and 
two  years  are  necessary  to  produce  a  new  spacecraft  and  to  procure  a  new 
launch  vehicle  respectively.  For  the  1  sat  / 1  set  configuration,  there  is  a  2.9% 
probability  of  three  failures  occurring  within  one  year  of  each  other.  Thus  three 
spares  and  three  sets  of  instruments  are  required  at  the  beginning  of  the 
mission  to  ensure  achievement  of  mission  availability  requirements  with  99% 
confidence.  Adding  redundant  instruments  to  the  single  satellite  lowers  this 
probability  of  three  failures  within  the  year  to  1.1%.  Therefore  the  1  sat  /2  set 
configuration  also  requires  three  spare  spacecraft  at  the  beginning  of  the 
mission. 
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Only  three  launch  vehicle  spares  are  required  for  the  1  sat  /2  set  configuration 
since  the  probability  of  four  failures  within  the  two  years  necessary  to  procure  a 
new  launch  vehicle  drops  below  1%  with  the  addition  of  redundant  sensors. 
The  cost  savings  from  needing  one  less  spare  launch  vehicle  is  not  enough  to 
offset  the  additional  cost  of  the  redundant  sensors. 

With  a  full  on-orbit  complement  of  nine  smaller  satellites,  the  3  sat  /  1  set 
configuration  requires  only  four  spare  buses  and  three  spare  sets  of 
instruments.  This  results  in  a  large  decrease  in  the  cost  of  the  spacecraft  spares 
as  compared  to  the  single  satellites  configurations.  Because  there  are  nine 
satellites,  seven  launch  vehicle  spares  must  be  available  to  ensure  confidence 
of  achieving  mission  availability.  These  launch  vehicles,  LMLV  lls,  cost  $21 
million  compared  with  the  $60  million  cost  to  procure  spare  Delta  lls  for  the 
single  satellite  configuration.  Overall  ,  the  cost  to  provide  necessary  spares  is 
lowest  for  the  3  sat  / 1  set  configuration. 

Distributing  the  primary  instruments  among  three  satellites  significantly 
increases  the  reliability  of  each  individual  satellite.  While  this  effect  was 
apparent  in  the  number  of  spares  required  it  also  influences  the  cost  of 
expected  replacements  procured  during  the  mission  life.  Higher  satellite 
reliability  and  lower  launch  costs  to  replace  the  satellites  that  do  fail  results  in 
the  3  sat  /  1  set  configuration  having  the  lowest  expected  replacement  cost. 
Once  again  the  slight  increase  in  reliability  gained  from  adding  redundant 
primary  instruments  for  the  3  sat  /  2  set  configuration  is  outweighed  by  the 
higher  bus,  payload,  and,  launch  costs. 

The  combined  effect  of  these  three  categories  of  cost  is  shown  by  the  total  ten 
year  non-operations  costs  data  in  Figure  2-3.  The  3  sat  /  1  set  configuration 
appears  less  expensive  than  the  1  sat  / 1  set  configuration. 
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Figure  2-4  displays  the  same  data  as  in  Figure  2-3  although  broken  down  into 
RDT&E,  payload,  bus,  and  launch  costs,  expended  during  the  10  year  mission. 
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Figure  2-4  Breakdown  of  costs  over  10  year  mission  life 


2.4.3  Sensitivity  of  Costs  to  Instrument  Reliability 

The  sensitivity  of  the  10  year  mission  cost  to  sensor  reliability  also  indicates  the 
best  choice.  The  previous  analysis  assumed  a  sensor  reliability  of  0.86  over  7 
years.  If,  however,  the  reliability  of  the  sensors  were  to  fall  short  of  the  stated 
goal;  for  example  0.70  over  seven  years,  the  10  year  mission  cost  for  the  3  sat  / 
1  set  and  1  sat  / 1  set  configurations  would  be  as  shown  in  Figure  2-5.  Because 
the  primary  instruments  are  distributed  across  three  satellites,  a  decrease  in 
their  reliability  results  in  only  a  relatively  small  increase  in  cost  for  the  3  sat  / 1 
set  configuration.  Because  all  of  the  primary  sensors  are  on  a  single  satellite  in 
the  1  sat  / 1  set  configuration,  a  decrease  in  sensor  reliability  significantly 
reduces  the  overall  reliability  of  that  satellite  leading  to  a  jump  of  over  $500 
million  in  the  10  year  mission  cost.  The  more  prudent  choice  to  hedge  against 
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the  risk  of  a  large  cost  increase  due  to  lower  than  expected  sensor  reliability  is 
the  3  sat  / 1  set  configuration. 


Figure  2-5  10  year  mission  cost  versus  sensor  reliability 


Distribution  of  sensors  among  several  satellites  is  appropriate  when  the 
savings  from  increased  reliability  outweigh  the  increased  cost  to  deploy  the 
distributed  system.  For  satellites  designed  to  a  short  life  (e.g.  3  years),  it  is  more 
likely  that  the  bus  will  reach  the  end  of  its  propellant  supply  or  battery  life  than 
for  a  sensor  to  fail.  Thus,  any  gains  in  reliability  from  distribution  are  irrelevant 
and  the  distributed  system  costs  more  over  the  mission  life.  Further,  distribution 
is  less  advantageous  for  systems  with  high  sensor  reliability.  As  sensor 
reliability  increases,  the  overall  reliability  of  the  1  sat  / 1  set  configuration  quickly 
approaches  that  of  the  3  sat  / 1  set  configuration.  The  high  initial  deployment 
costs  of  the  3  sat  /  1  set  configuration  can  no  longer  be  justified  when 
replacement  costs  for  the  1  sat  / 1  set  configuration  are  on  par  with  those  of  the 
3  sat  /  1  set  configuration. 
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It  should  be  mentioned  that  no  attempt  was  made  to  account  for  a  decrease  in 
the  procurement  cost  for  lower  reliability  sensors.  Figure  2-5  indicates  that  only 
a  small  increase  in  the  3  sat  / 1  set  mission  cost  takes  place  when  sensors 
having  a  reliability  of  0.7  over  7  years  are  used.  Had  the  decrease  in  sensor 
procurement  cost  also  been  factored  in,  this  small  increase  may  have  been 
further  reduced. 

Table  2-1 1  summarizes  the  four  candidate  cluster  configurations’  influence  on 
the  measures  of  effectiveness.  Each  of  the  four  configurations  provide  an 
opportunity  for  the  system  to  measure  required  environmental  data  with 
necessary  co-registration  requirements.  Because  increased  schedule  risk  is 
due  primarily  to  the  use  of  unproven  technologies,  the  choice  of  how  the 
sensors  are  distributed  neither  significantly  increases  nor  decreases  the  risk  of 
the  schedule  slipping.  The  configurations  with  the  redundant  set  of  instruments 
were  shown  to  greatly  increase  cost  while  minimum  cost  was  achieved  by  the 
configurations  with  only  a  single  set  of  instruments.  Only  the  3  sats  /  1  set 
configuration,  however,  protects  against  a  large  increase  in  cost  due  to  lower 
than  expected  sensor  reliability. 

Table  2-11  Comparison  with  respect  to  Measures  of  Effectiveness 


Weight 

Factor 

1  Sat 

1  Set 

1  Sat 

2  Sets 

3  Sats 

1  Set 

Minimize  Cost 

4 

X 

X 

✓ 

Minimize  Cost  Risk 

3 

X 

S 

✓ 

Minimize  Schedule  Risk 

1 

Q 

n 

o 

Provide  Opportunity  To 

Meet  Schedule  Reqs. 

8 

✓ 

✓ 

✓ 

■  -  Neutral  S  -  Allows  M 

IOE  *  -  Prevents  MOE 

2.5  OTHER  IMPLEMENTATION  CONSIDERATIONS 

Current  levels  of  automation  in  satellite  systems  reflect  an  incremental  evolution 
that  is  based  on  a  high  level  of  human  involvement.  Historically,  this  has  been 
a  result  of  the  desire  to  reduce  risk  and  due  to  limited  technological  capabilities. 
Due  to  the  dependence  on  humans  to  perform  tasks,  operations  costs  can  make 
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up  a  significant  portion  of  the  life  cycle  costs  for  a  satellite  system  [8].  In 
addition,  human  error  continues  to  be  a  major  cause  of  spacecraft  anomalies 
and  failures.  With  the  introduction  of  large  constellations  or  clusters  of  satellites, 
some  automation  of  operations  will  be  required  to  reduce  costs  while 
maintaining  availability  (the  probability  of  meeting  system  requirements  at  a 
given  time)  [9],  [10],  [11].  Despite  recognition  of  the  need,  there  is  reluctance  to 
automate.  The  large  investments  and  high  risks  involved  in  space  ventures  has 
lead  to  a  conservative  industry.  In  addition,  the  desire  to  reduce  cycle  times  for 
new  programs  also  favors  significant  re-use  of  proven  technologies  and  thus, 
low  levels  of  automation. 

Figure  2-6  qualitatively  represents  the  cost  and  availability  characteristics  of  a 
hypothetical  satellite  system  with  respect  to  an  increasing  level  of  automation. 
As  low  levels  of  automation  are  introduced  into  the  system,  the  operating  costs 
decrease,  principally  due  to  a  decrease  in  the  number  of  human  operators 
(Figure  2-6A).  At  some  point,  however,  the  increases  in  design  and 
development  costs  due  to  software  development  outweigh  the  decrease  in 
operation  costs. 

As  shown  in  Figure  2-6B,  availability  may  decrease,  increase,  or  be  unaffected 
by  the  level  of  automation.  For  tasks  that  are  simple,  well  understood,  or 
periodic  (such  as  routine  station-keeping  on  a  geostationary  satellite), 
availability  may  increase  with  increasing  automation  (Task  A).  This  is  true  in  the 
cases  when  human  errors  are  more  likely  than  software  errors,  or  when  the 
impact  of  unanticipated  situations  is  negligible.  For  complex,  rare,  or 
unexpected  functions,  availability  may  decrease  as  humans  are  removed  from 
the  loop  (T ask  C).  This  can  occur  when  the  automation  is  unable  to  resolve 
problems  that  could  have  been  resolved  by  a  human  or  when  the  automation 
fails  to  accurately  inform  the  human  of  the  situation.  There  may  be  some 
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functions  for  which  the  availability  is  nearly  independent  of  the  level  of 
automation  (Task  B). 

Development  and  Availability  Life  Cycle  Costs 


A  B  C 

Figure  2-6:  Effect  of  Automation  on  Cost  and  Availability 

An  increase  in  system  availability  translates  into  an  increase  in  revenues  for  a 
commercial  system,  or  an  increased  ability  to  perform  an  objective  for  science 
or  military  systems.  Thus,  the  potential  for  failure  can  be  represented  by  an 
opportunity  cost  which  represents  revenues  forgone  as  a  result  of  increased 
system  down-time.  For  commercial  systems,  the  opportunity  cost  can  be  added 
to  the  development  and  operations  costs  to  form  the  life  cycle  cost. 
Determination  of  opportunity  costs  requires  additional  data  such  as  the 
relationship  between  a  particular  function  and  revenue.  For  science  or  military 
applications,  the  definition  of  an  opportunity  cost  may  be  difficult.  In  such  cases, 
the  development  and  operations  cost  would  be  compared  against  the 
availability  without  attempting  to  define  the  life  cycle  cost.  The  combination  of 
the  two  curves  Figure  2-6A  and  Figure  2-6B  are  represented  in  Figure  2-6C, 
showing  the  overall  life  cycle  cost  which  is  defined  as  the  sum  of  development, 
operations,  and  opportunity  costs. 

In  the  example  of  Figure  2-6,  there  exists  an  optimum  level  of  automation  at 
which  life  cycle  cost  is  minimized.  Due  to  the  complexity  of  the  satellite  system, 
a  methodology  is  needed  that  can  model  the  effect  that  automation  has  on  costs 
and  system  availability.  Such  tools  would  enable  system  engineers  to  identify 
those  functions  that  should  be  automated. 
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Constellations  of  satellites  introduce  a  twist  in  the  automation  trade  since 
automation  is  usually  associated  with  high  non-recurring  costs.  As  the  level  of 
automation  increases,  operations  costs  are  expected  to  decrease,  and 
development  costs  are  expected  to  increase.  As  the  number  of  satellites 
increases  for  a  given  level  of  automation,  the  operations  costs  will  increase 
linearly.  Development  costs,  however,  will  lag  a  linear  relationship  due  in  part 
to  a  learning  curve,  and  in  part  to  the  much  lower  recurring  costs  of  automation 
relative  to  the  non-recurring  cost.  Thus,  functions  which  are  not  suitable  for 
automation  for  a  single  satellite  may  be  desirable  for  a  constellation.  In  fact, 
crossover  points  can  be  identified  which  define  constellation  sizes  over  which 
certain  levels  of  automation  are  optimal. 

Constellations  of  microsatellites  will  likely  contain  a  very  large  number  of 
satellites  in  order  to  perform  a  useful  mission.  As  discussed  above,  for  a  given 
level  of  automation,  as  the  constellation  size  increases,  operations  costs  will 
scale  linearly  while  development  costs  will  take  advantage  of  economies  of 
scale.  Therefore,  the  operations  costs  for  these  large  constellations  will 
dominate  the  life  cycle  costs,  and  high  levels  of  automation  will  be  required  for 
many  functions. 


2.6  SUMMARY  AND  CONCLUSIONS 

A  distributed  architecture  makes  sense  if  it  can  offer  reduced  cost  or  improved 
performance.  Functional  requirements  specify  minimum  levels  of  acceptable 
performance,  and  include  resolution,  rate,  integrity  and  availability 
requirements.  Viable  systems  must  satisfy  these  requirements  throughout  their 
lifetime.  Compensation  must  be  made  following  failures  that  cause  a  violation  of 
requirements.  “Improved  Performance”  thus  relates  to  the  ability  of  the  system 
to  satisfy  requirements  with  a  higher  probability.  “Reduced  Cost”  corresponds 
to  lower  lifetime  costs  that  include  the  expected  failure  compensation  costs. 
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Because  the  performance  requirements,  and  the  associated  probability  of 
satisfying  them,  are  embedded  in  its  calculation,  lifetime  cost  is  a  useful  metric 
for  architecture  analysis. 


Distribution  can  offer  improvements  in  isolation  (resolution),  rate,  integrity  and 
availability.  The  improvements  are  not  all-encompassing,  and  in  many  cases 
are  application  specific.  Nevertheless,  it  appears  that  adopting  a  distributed 
architecture  can  result  in  substantial  gains  compared  to  traditional  deployments. 
Some  of  the  more  important  advantages  that  distribution  may  offer  are: 

•  Improved  resolution  corresponding  to  the  large  baselines  that  are  possible 
with  widely  separated  antennas  on  separate  spacecraft  within  a  cluster. 

•  Higher  net  rate  of  information  transfer,  achieved  by  combining  the  capacities 
of  several  satellites  in  order  to  satisfy  the  local  and  global  demand. 

•  Improved  availability  through  redundancy  and  path  diversity.  Frequently,  the 
cost  of  adding  a  given  level  of  redundancy  is  less  for  a  distributed 
architecture. 

•  Lower  failure  compensation  costs  due  to  the  separation  of  important  system 
components  among  many  satellites;  only  those  components  that  break  need 
replacement. 

There  are  some  problems,  specific  to  distributed  systems  of  small  satellites,  that 
must  be  solved  before  the  potential  of  distributed  architectures  can  be  fully 
exploited.  The  most  notable  of  these  problems  are: 

•  An  increase  of  system  complexity,  leading  to  long  development  time  and 
high  costs 

•  Inadequacy  of  the  data  storage  capacity  that  can  be  supported  by  the 
modest  small  satellite  bus  resources 

•  Difficulty  of  maintaining  signal  coherence  among  the  apertures  of  separated 
spacecraft  arrays,  especially  when  the  resolution  requirements  are  high,  or 
the  target  is  highly-dynamic. 

The  resolution  of  these  issues,  and  the  proliferation  of  microtechnology,  could 
lead  toward  a  drastic  change  in  the  satellite  industry.  It  seems  clear  that 
distribution  offers  a  viable  and  attractive  alternative  for  some  missions.  Large 
constellations  of  hundreds  or  thousands  of  small-  and  micro-satellites  could 
feasibly  perform  almost  all  of  the  missions  currently  being  carried  out  by 
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traditional  satellites.  For  some  of  those  missions,  the  utility  and  suitability  of 
distributed  systems  looks  very  promising.  More  analysis  is  warranted  in  order  to 
completely  answer  the  question  of  where  and  when  distribution  is  best  applied, 
but  the  potential  prospects  of  huge  cost  savings  and  improvements  in 
performance  are  impossible  to  ignore.  It  therefore  seems  inevitable  that 
massively  distributed  satellite  systems  will  be  developed  in  both  the  commercial 
and  military  sectors.  We  are  living  in  a  time  of  great  changes,  and  the  space 
industry  has  not  escaped.  Over  the  last  few  years,  “faster,  cheaper,  better”  has 
been  the  battle  cry  of  those  engineers  and  administrators  trying  to  instigate 
changes  to  improve  the  industry.  “Smaller,  modular,  distributed”  may  be  their 
next  verse. 


References 

1.  “Foresight”  Final  Report  from  MIT  16.89  Space  Systems  Engineering 
class,  Massachusetts  Institute  of  Technology,  Cambridge,  MA,  May,  1997. 

2.  Shaw,  G.,  Yashko,  G.,  Schwarz,  R.,  Wickert,  D.  and  Hastings,  D.  “Analysis 
Tools  and  Architecture  Issues  for  Distributed  Satellite  Systems”  in 
Microenqineerinq  for  Aerospace  Systems  H.  Helvijan  (ed),  Aerospace 
Press,  El  Segundo,  CA  (1997). 

3.  Infotech  Report,  System  Reliability  and  Integrity,  Infotech  International 
Limited,  1978 

4.  Request  For  Proposal  #F04701  -95-R-0032.  NPOESS  program  office,  Feb. 
1996. 

5.  Larson,  Wiley  J.  and  Wertz,  James  R.  Space  Mission  Analysis  and  Design 
2nd  ed.,  Microcosm,  Inc.  and  Kluwer  Academic,  Torrance,  California,  1992 

6.  Larson  and  Wertz.  Space  Mission  Analysis  and  Design.  Microcosm  Press, 
1994. 

7.  Greenberg,  Joel.  AIAA  Paper  96-1 1 12-CP,  1996. 

8.  Congor,  Robert  Microcosm  Autonomous  Navigation  System  (MANS). 
Presentation,  August,  1995. 


44 


9.  Tandler,  John.  “Automating  the  Operations  of  the  ORBCOMM 
Constellation”,  10th  Annual  AIAA/USU  Conference  on  Small  Satellites, 
Utah,  September  1996 

10.  Hornstein,  Rhoda  Shalter.  “On-Board  Autonomous  Systems:  Cost 
Remedy  for  Small  Satellites  or  Sacred  Cow?”  46th  International 
Astronautical  Congress,  Oslo,  Norway,  Oct.  2-6,  1995. 

11.  Collins,  John  T,  Simon  Dawson,  &  James  R.  Wertz.  “Autonomous 
Constellation  Maintenance  System.”  10th  Annual  AIAA/USU  Conference 
on  Small  Satellites,  1996. 


45 


This  page  intentionally  left  blank 


46 


Chapter  3 

Thruster  Requirements  for  Local 
Satellite  Clusters 


This  chapter  examines  the  propulsive  requirements  necessary  to  maintain  the  relative 
positions  of  satellites  orbiting  in  a  local  cluster.  Formation  of  these  large  baseline 
arrays  could  allow  high  resolution  imaging  of  terrestrial  or  astronomical  targets  using 
techniques  similar  to  those  used  for  decades  in  radio  interferometry.  A  key  factor  in 
the  image  quality  is  the  relative  positions  of  the  individual  apertures  in  the  sparse 
array.  The  relative  positions  of  satellites  in  a  cluster  are  altered  by  “tidal”  accelerations 
which  are  a  function  of  the  cluster  baseline  and  orbit  altitude.  These  accelerations 
must  be  counteracted  by  continuous  thrusting  to  maintain  the  relative  positions  of  the 
satellites. 

3.1  INTRODUCTION  TO  CLUSTER  MISSIONS 

Interferometry  has  been  used  for  decades  to  produce  images  of  astronomical  objects 
in  radio  wavelengths  with  resolutions  rivaling  that  of  ground-based  optical  systems  [1]. 
This  technique,  when  deployed  across  a  cluster  of  satellites,  may  be  used  to  produce 
high  resolution  remote  sensing  images  of  terrestrial  targets  from  space  or  provide  a 
platform  above  the  atmosphere  for  space-based  viewing  of  astronomical  objects. 


3.1.1  Synthetic  Apertures 

The  diffraction  limited  ground  resolution  of  typical  filled  apertures  is  given  by  [2] 
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(3-1) 


^  ~ 


b 


where  b  =  typical  aperture  dimension  and  Rs=  slant  range  to  target. 


Equation  (3-1)  indicates  that  increasing  the  aperture  size  improves  the  resolving 
power  of  the  instrument.  Aperture  size  is  limited,  however,  by  size  and  weight 
constraints  of  the  launch  vehicles  which  place  the  satellites  in  orbit.  Synthetic 
Aperture  Radar  (SAR)  can  improve  resolution  at  RF  wavelengths  from  Low  Earth  Orbit 
(LEO)  by  utilizing  the  Doppler  shifting  of  signals  to  “synthesize”  an  aperture  much 
larger  than  the  satellite’s  filled  aperture.  Resolutions  of  SAR  systems  are  given  by 

(3-2) 

where  b  is  again  a  characteristic  dimension  of  the  filled  aperture.  Contrary  to  intuition, 
equation  (3-2)  indicates  that  reducing  the  filled  aperture  size  improves  the  resolving 
power  of  the  radar.  At  some  point,  however,  aperture  temperature,  power  output,  and 
gain  limitations  preclude  the  use  of  smaller  apertures.  Table  3-1  lists  allowable 
resolutions  for  a  typical  10  m  filled  aperture  sensing  in  the  optical,  infrared  (IR),  or 
radio  frequency  (RF)  wavelengths  at  several  altitudes. 


able  3-1  Current  Achievable  Resolution  with  10m  Aperture 


Altitude 

Optical 

X=0.5pm 

Infrared  (IR) 
X=10  pm 

Radio  (RF) 
k=3  cm 

1 ,000  km 

5  cm 

1  m 

5  m  (SAR) 

1 0,000  km 

50  cm 

10  m 

30  km  1 

35,768  km 

1.8  m 

18  m 

107  km 

3.1.2  Interferometry  /  Sparse  Arrays 

Another  method  to  increase  resolution  is  interferometry.  Interferometers  use  separate 
apertures  spaced  some  distance  apart  to  create  a  sparse  aperture  as  sketched  in 
Figure  3-1 
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Figure  3-1  Ground-based  Radio  Interferometer 


Interferometers  were  originally  used  in  radio  astronomy  to  distinguish  discrete  radio 
sources  from  the  diffuse  background.  Synthesis  is  obtained  by  observing  separately 
all  the  interferometer  pairs  that  exist  within  a  large  aperture.  A  chessboard  with  only 
two  chess  pieces  may  be  used  as  an  analogy.  All  combinations  of  the  pieces  on  the 
board  must  be  sampled  to  recreate  a  uniform  aperture.  To  accomplish  this,  the 
receiving  elements  are  designed  to  be  mobile  across  the  so-called  u-v  plane.  As 
baselines  increased,  the  rotation  of  the  earth  was  used  to  sweep  a  range  of  spacings 
with  fixed  antennae.  Very  large  baselines  composed  of  widely  separated  antennae 
soon  made  sampling  the  entire  u-v  plane  impractical.  Computer  calibration 
techniques  are  now  commonly  used  for  sparse  aperture  synthesis. 

The  position  and  spacing  of  the  elements  is  key  to  the  quality  of  the  image  produced. 
The  paths  traced  out  by  the  electro-magnetic  waves  must  be  carefully  controlled  so 
that  the  signals  may  be  coherently  combined.  This  tolerance  is  typically  A/20.  For  a 

regularly  spaced  array,  a  positional  error  of  a  few  wavelengths  can  cause  significant 
sidelobes,  even  if  the  faulty  positions  are  known  exactly  so  that  each  element  can  be 
correctly  phased  [3].  In  cases  where  the  relative  positions  of  the  element  are  not 
known  (such  as  with  random  arrays),  calibrations  can  be  carried  out  to  improve  the 
quality  of  the  images  providing  the  relative  positions  do  not  change  from  one  image  to 
the  next.  Although  routine  at  RF  wavelengths,  sparse  aperture  formation  at  IR  and 
visible  wavelengths  is  still  experimental. 


The  achievable  angular  resolution  is  diffraction  limited  in  accordance  with  the 
Rayleigh  criterion  (Eq.  3-1).  The  size  of  the  aperture,  however,  is  the  maximum  linear 
distance  between  the  individual  elements.  This  distance  is  known  as  the  baseline,  B, 
so  that  [4] 


Table  3-2  lists  the  angular  resolution  versus  sampled  wavelengths  for  four  existing 
interferometers  along  with  the  baseline  of  each  system.  Although  aperture  synthesis  is 
well  established  in  the  radio  spectrum,  imaging  of  astronomical  sources  with  multi¬ 
aperture  interferometers  is  just  beginning  with  the  construction  of  the  Infrared-Optical 
Telescope  Array  (IOTA)  at  the  F.L.  Whipple  Observatory  on  Mount  Hopkins,  AZ.  [5] 


Table  3-2  Angular  resolution  of  existing  interferometers 


Interferometer 

Baseline 

X 

Ang.  Resol.  | 

VLBI 

8,500  km 

Icm-lm 

1  e-4”  to  1  e-2” 

VLA 

35  km 

1cm-5m 

0.1”  to  20” 

Westerbrook 

3.2km 

lOcm-lm 

5”  to  100” 

Merlin 

135  km 

lOcm-lm 

0.1”  to  1”  1 

Source:  Wohllenben  [6] ,  Pg  17 


3.1.3  Local  Satellite  Clusters 

Forming  space-based  sparse  apertures  could  be  accomplished  through  the  use  of  a 
local  satellite  cluster.  Satellites  in  a  local  cluster  orbit  in  formations  such  as  those 
shown  in  Figure  3-2  and  Figure  3-3.  Although  one  or  more  reference  satellites  will  be 
in  standard  Keplerian  (ie.  inertial)  orbits,  maintaining  the  formation  will  require  the 
other  satellites  to  orbit  in  planes  parallel  to  the  reference  orbits.  These  non-inertial 
orbits  are  characterized  by  either  a  focus  which  is  not  located  at  the  Earth’s  center  of 
mass  (Figure  3-2),  or  orbital  velocities  which  do  not  provide  the  proper  centripetal 
acceleration  to  offset  gravity  at  that  altitude  (Figure  3-3).  As  expected,  the  Earth’s 
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gravitation  will  act  to  move  these  satellites  into  Keplerian  orbits.  It  will  be  shown  that 
continuous  low-level  thrusting  is  required  to  maintain  each  satellite’s  position. 


Figure  3-2  Cluster  Formation  Normal  to  Reference  Orbit  Plane 


non-inertial  orbit 


Figure  3-3  Cluster  Formation  in  Plane  of  Reference  Orbit 

By  coherently  adding  the  signals  received  by  several  satellites,  the  cluster  would,  in 
effect,  create  a  sparse  aperture  many  times  the  size  of  a  real  aperture.  The  sparse 
aperture’s  large  size  could  vastly  increase  the  level  of  resolution  possible.  There  are 
several  possible  advantages  of  creating  a  sparse  array  cluster.  For  example: 

1)  Extensive  earth  coverage  could  be  achieved  at  GEO  with  a  resolution  similar  to 
that  of  current  satellites  in  LEO. 
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2)  By  turning  satellites  within  the  cluster  “on  and  off”  it  may  be  possible  to  alter  the 
dimensions  of  the  sparse  aperture  and  thus  zoom  in  on  a  target  detected  in  a 
broad,  coarse  field  of  view. 

The  gravitational  forces  which  act  on  the  satellites  in  non-inertial  orbits  will  now  be 
examined  in  more  detail. 


3.2  ORBITAL  DYNAMICS  OF  SATELLITE  CLUSTERS 
3.2.1  Tidal  Forces 

Figure  3-4  shows  the  coordinate  system  for  a  simple  two  satellite  cluster.  Satellite  1  is 
in  an  inertial  reference  orbit  at  an  distance,  R,  from  the  center  of  the  Earth.  Satellite  2’s 
position  relative  to  satellite  1  is  r. 


Figure  3-4  Satellite  Cluster  Coordinate  System. 


The  motion  satellite  2  with  respect  to  1  can  be  found  using  the  linearized  equations  of 
relative  motion  for  a  circular  reference  orbit.  These  equations,  known  as  the 
Clohessey-Wiltshire  equations  [7],  are 
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where  Q  =  ^p/R3. 


One  would  intuitively  expect  a  rz  term  due  to  displacements  along  y.  It  is  eliminated 
by  the  linearization  process  in  the  derivation  of  equation  (3-4)  although  it  does  appear 
in  the  full  length  equations. 


r 

r  = 


v 


z-2Qx) 


'  n2x 
n2y 

S22(R+z) 


(3-5) 


where  c= 


x2  +  y2  +  z2>\ 

*2  y 


Section  3  identifies  the  maximum  displacements  from  the  reference  orbit  that  are 
feasible  due  spacecraft  mass,  power,  and  volume  constraints.  For  those  size  clusters, 
this  extra  Tz  term  is  negligible. 

Equation  (3-4)  when  multiplied  by  the  m^,  describes  the  required  on-board  thrusting 
necessary  to  produce  desired  accelerations  and  velocities  relative  to  some  reference 
satellite  in  an  inertial  orbit.  Notice  that  for  x  =  y  =  z  =  x  =  y  =  z  =  x  =  y  =  z  =  0,  (3-4) 
reduces  to  the  case  of  a  satellite  in  a  inertial  orbit  where  no  thrusting  is  necessary  to 
maintain  its  orbit  under  two-body  orbit  assumptions. 


Even  with  x  =  y  =  z  =  x  =  y  =  z  =  0,  y  and  z  displacements  create  “tidal”  forces  which 
will  tend  to  move  the  satellite  into  an  inertial  orbit  as  described  by  Janson  [8].  These 
tidal  forces  need  to  be  counteracted  by  thrusting  in  order  to  maintain  the  cluster.  Tidal 
forces  arise  from  displacements  along  z  because  the  cluster  satellite  is  constrained  to 
orbit  with  the  same  velocity  as  the  reference  satellite  yet  at  a  different  altitude.  Thus, 


unlike  the  reference  satellite,  the  gravitational  attraction  of  the  Earth  is  not  exactly 
offset  by  the  centripetal  acceleration  due  to  the  satellite’s  circular  motion. 
Displacements  along  y  force  the  satellite  to  orbit  in  a  plane  parallel  to  reference 
satellite.  In  that  case,  only  a  component  of  the  gravity  vector  lies  in  the  plane  of  the 
orbit.  Once  again,  the  centripetal  acceleration  is  not  offset  by  Earth’s  gravity.  A 
displacement  along  x  can  be  considered  a  displacement  along  z  occurring  at  a  point 
further  ahead  in  the  reference  orbit. 


3,2.2  Impulsive  vs.  Continuous  Thrusting 

To  remain  within  the  allowable  relative  position  tolerances  discussed  previously,  the 
cluster  of  satellites  could  use  either  a  series  of  impulsive  thruster  firings  at  regular 
intervals  or  apply  a  lower,  continuous  thrust. 

Consider  the  case  of  a  satellite  displaced  along  y  as  shown  in  Figure  3-5.  If  no  thrust 
were  applied  in  the  y  direction,  the  displaced  satellite  would  oscillate  around  the 
reference  orbit  completing  one  oscillation  every  orbit.  In  essence,  it  would  be  in  an 
inertial  orbit  inclined  with  respect  to  the  reference  orbit.  This  motion  relative  to  the 
reference  orbit  would  temporarily  move  the  satellite  beyond  the  allowed  tolerance  for 
elements  of  a  sparse  array.  Therefore,  the  thruster  will  need  to  fire  an  impulse  bit 
which  keeps  the  satellite  within  the  tolerance  region. 
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Figure  3-5  Impulsive  Thrusting  Procedure 


At  time  t=0,  halfway  between  impulse  firings,  x  =  z  =  x  =  z  =  x  =  z  =  rx=ry=rz=0. 
Equation  (3-4)  then  reduces  to 

y  +  Q2y  =  0  (3-6) 

This  undamped  linear  oscillator  has  the  solution 

y(t)  =  y0  cos(Qt)  (3-7) 


y(0  =  -y0^sin(Qt) 


(3-8) 


y(t)  =  -yon2cos(Qt) 


(3-9) 


Then,  from  equation  (3-7)  the  time  at  which  the  satellite  would  drift  beyond  the 
tolerance  limit  is 


1 

t  =  — cos 
Q 


-1 


Vo -top 
<  y0  ; 


(3-10) 


Figure  3-6  shows  the  plot  of  equation  (3-10)  at  several  altitudes.  It  can  be  seen  from 
the  figure  that,  at  all  altitudes,  as  the  required  relative  position  tolerance  decreases  the 
time  the  satellite  is  able  to  drift  before  firing  an  impulse  decreases.  The  largest  value 
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of  tolerance/displacement  will  be  for  RF  sensing  sparse  apertures  in  LEO.  For  an  RF 
cluster  satellite  sensing  at  k=.05m,  tolerances  will  be  on  the  order  of  2.5  mm.  At  1 000 

km  altitude  and  displaced  25m,  this  tolerance/displacement  is  1e-4.  Figure  3-6 
indicates  that  the  satellite  could  drift  for  14  seconds.  In  the  optical  region,  however, 
allowable  drift  time  is  approximately  0.05  seconds. 


Figure  3-6  Time  Until  Tolerance  is  Exceeded 


Assuming  that  the  satellite’s  y  velocity  after  the  impulse  is  exactly  negative  that  before 
the  impulse,  and  that  the  tidal  forces  and  impulse  thrust  are  constant  during  the  time 
the  impulse  is  applied,  then 

=  2y0Qsin(Qt)  +  [y0Q2  cos(Qt)j(At)  (3-1 1 ) 

ms/ c 

where  t  is  found  from  equation  (3-10)  and  At  is  the  thruster  firing  time. 


The  total  impulsive  firings  during  the  mission  life  is  then 


Nf  = 


life 

2t  + At 


(3-12) 
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Calculating  the  total  impulsive  AV  expended  over  the  mission  life, 


AV  =  Nf 

ms/c 


(3-13) 


Figure  3-7  compares  the  AV  expended  for  continuous  thrusting  to  that  of  impulse 
thrusting.  The  figure  shows  that  for  a  given  displacement,  as  the  tolerance  region  in 
which  the  satellite  can  drift  increases,  only  negligible  amounts  of  AV  can  be  saved  by 
using  impulsive  thrusts.  As  the  tolerances  tighten  (and  the  time  between  impulsive 
thrusts  shortens),  the  AV  expended  for  a  impulsive  thrust  approaches  that  expended 
by  continuous  thrusting. 
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Figure  3-7  Propellant  Savings  from  Impulsive  Thrusting. 


The  previous  analysis  would  indicated  that,  at  the  broadest  relative  position  tolerance 
levels,  the  maximum  time  between  required  impulsive  thrusts  is  on  the  order  of 
seconds.  For  optical  and  infrared  wavelengths,  this  time  is  on  the  order  of  tenths  of  a 
second.  Further,  at  these  short  thrusting  intervals,  only  minimal  amounts  of  propellant 
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savings  result  from  impulsive  thrusting.  Therefore,  continuous  thrusting  will  be 
assumed  for  the  rest  of  this  analysis. 


3.2.3  Dynamic  Clusters 

For  the  situation  where  satellite  positions  are  fixed  relative  to  each  other,  the  total  AV 
expended  during  the  satellite’s  mission  life  is 

AV  =  (f)life  (3-14) 

4- 

where  T  is  the  acceleration  produced  continuously  by  an  on-board  propulsive  system 
to  counteract  tidal  forces.  Equation  (3-4)  indicates  that  the  required  thrust  per  unit  of 
spacecraft  mass  is  larger  for  a  satellite  stationed  at  a  greater  distance  from  the 
reference  orbit.  From  equation  (3-14)  for  a  given  mission  life,  a  satellite  stationed  far 
from  the  reference  orbit  will  consume  more  fuel  than  a  satellite  stationed  closer.  It  may 
be  desirable,  therefore,  to  rotate  the  positions  of  the  satellites  during  their  life  as  a 
means  of  distributing  the  fuel  consumption  amongst  all  the  satellites  in  the  cluster. 

Consider  the  case,  illustrated  in  Figure  3-8,  where  the  satellites  in  a  cluster  continually 
rotate  in  a  circular  pattern  relative  to  the  reference  satellite.  Figure  3-8  shows  only  two 
of  the  many  satellites  comprising  the  cluster.  This  is  one  of  many  possible  methods  by 
which  satellites  in  a  cluster  change  position,  but  it  serves  well  as  an  example  of  a 
“dynamic  cluster”. 


Figure  3-8  Dynamic  Cluster  Formation 
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For  this  situation, 


x  =  rsin(0) 
y  =  rcos(0) 
z  =  0 


(3-15) 


By  substituting  the  relations  from  equation  (3-15),  equation  (3-4)  becomes 

rx  =  -i02sin0 

ry  =  -rcosO^Q2  +  02  j  (3-1 6) 

Tz  =  2r(cos0)Q0 

where  Q  =  -y/jx/R3 . 

It  can  be  seen  from  equation  (3-16)  that  the  required  thrust  levels  of  a  satellite  in  a 
dynamic  cluster  are  a  function  of  the  cluster  radius,  r,  the  cluster’s  reference  altitude, 
the  position  of  the  satellite,  0,  as  well  as  the  rate  at  which  the  satellite  rotates  about  the 

reference  satellite,  0.  Figure  3-9  shows  AV  expended  by  all  satellites  in  the  example 

dynamic  cluster  can  be  as  low  as  54%  that  expended  by  a  satellite  at  the  far  edge  of  a 
static,  non-rotating  cluster. 


Number  of  Rotations  During  Mission  Life 

Figure  3-9  Reduction  in  Cluster  Maintenance  AV 
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As  the  number  of  rotations  about  the  reference  satellite  increases  beyond  about  1/3, 
the  AV  ratio  begins  to  increase  again  because  the  satellite  is  spending  more  time 
displaced  further  from  the  reference  orbit  along  -y.  At  a  rate  of  1  rotation  during  the 
mission,  the  AV  expended  is  about  65%  that  of  the  worst  case  satellite  in  a  static 
cluster.  The  AV  ratio  levels  off  at  this  value  for  more  than  5  rotations  during  the 
mission.  Although  not  shown  in  Figure  3-9,  as  the  rotation  rate  is  increased  further,  the 
AV  required  just  to  maintain  the  circular  motion  around  the  reference  satellite  begins  to 

dominate.  At  rate  of  a  few  hundred  rotations  during  the  mission,  the  AV  ratio  increases 
above  unity  indicating  it  is  no  longer  beneficial  to  rotate  the  cluster  satellites. 

It  should  be  stated  that  a  sparse  aperture  might  be  formed  using  a  constellation  of 
satellites  rather  than  a  local  cluster,  depending  on  whether  the  positional  tolerance 
problem  could  be  solved.  For  that  situation,  all  satellites  would  orbit  in  inertial  orbits 
and  the  sparse  aperture  would  be  formed  using  whatever  satellites  are  available  as 
they  fly  over  the  region  of  interest.  The  AV  for  maintaining  that  “sparse  aperture”  would 
be  zero.  The  analysis  throughout  the  rest  of  this  paper,  however,  is  focused  only  on 
the  requirements  of  maintaining  a  local  static  cluster. 


3.3  PROPULSION  SYSTEM  REQUIREMENTS 

Based  on  the  acceleration  levels  presented  in  section  2,  the  propulsion  system 
requirements  to  maintain  a  local  satellite  cluster  can  be  calculated.  The  feasible  range 
of  specific  impulse  and  efficiency,  constrained  spacecraft  mass,  volume  and  power 
considerations,  is  the  primary  characteristic  to  be  determined. 


3.3.1  Current  Thruster  Characteristics 

Two  classes  of  propulsion,  chemical  and  electric,  are  currently  used  as  station¬ 
keeping  thrusters  on-board  spacecraft.  These  two  classes  along  with  ranges  of 
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specific  impulse  and  efficiency  can  be  displayed  graphically  as  shown  in  Figure  3-10. 
These  ranges  reflect  some  technology  extrapolations  appropriate  for  the  next  5  to  1 0 
years.  For  exampe  PPT  efficiencies  are  currently  around  10%  but  work  is  under  way 
which  aims  to  extend  them  beyond  20%.  This  format  will  be  used  as  a  template  for  a 
comparison  of  current  thruster  capabilities  to  those  required  for  maintaining  cluster 
formations.  It  is  important  to  point  out  that  Figure  3-10  simply  shows  the  range  of  Isp 
vs.  ii.  If  a  local  cluster  formation  requires  a  combination  of  Isp  and  ii  which  falls  within 
the  shaded  box,  there  is  no  guarantee  a  thruster  with  those  characteristics  exists. 
However,  if  the  cluster  satellite  needs  a  combination  of  Isp  and  q  which  does  not  fall 

within  a  shaded  box,  then  it  can  be  said  that  no  current  thruster  exists  which  is  able  to 
meet  those  requirements. 


Efficiency 


Figure  3-10  Current  Thruster  Capabilities  [9],  [10] 

3.3.2  Methodology 

In  determining  the  feasible  range  of  thruster  Isp  and  ri  necessary  to  maintain  a  local 
satellite  cluster,  the  following  assumptions  were  made: 
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•  The  satellites  operate  as  a  static  cluster  maintaining  a  fixed  displacement  normal  to 
the  reference  orbit. 

•  Isp  and  r\  are  the  same  for  all  thrusters  used  in  cluster  maintenance  and  station¬ 
keeping 

•  Thrust  levels  of  each  thruster  are  constant  during  the  life  of  the  mission. 

•  Any  cluster  maintenance  and  station-keeping  maneuvers  are  performed 
continuously. 

•  The  spacecraft  and  propellant  tanks  are  spherical. 

•  Spacecraft  total  mass  is  constant  throughout  the  life  of  the  mission. 

•  An  additional  50m/s  of  propellant  per  year  is  allocated  for  traditional  station¬ 
keeping  requirements. 


A  spacecraft  displaced  a  given  distance  from  the  reference  orbit  using  a  thruster 
operating  at  a  given  t|  and  1^  must  satisfy  the  following  three  design  constraints. 

^tank  <  |  f-^tank 

In 

/max 


Ds/c 


mp+mpp 

ms/c 


Ds/c 

f 


mn  +m 


PP 


in 


^s/c 


V 

A 


ms/c 


(3-17) 


/max 


^s/c 


'max 


where 


mp  =  mass  of  propellant 
mpp=  mass  of  power  plant 


The  values  for  the  maximum  ratios  allow  margin  for  the  satellite  to  accomplish  tasks 
other  than  simply  operating  the  thrusters.  Specifically,  the  three  statements  gauge 
whether  or  not  there  is,  once  station-keeping  requirements  have  been  met,  sufficient 
volume,  mass,  and  power  available  on  board  the  satellite  to  perform  mission 
operations  such  as  payload,  communications,  etc.  If  the  answer  to  one  or  more  of 
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these  criteria  is  no,  then  it  is  concluded  that  a  satellite  of  mass,  m^,  and  propellant 
mass  fraction,  mp/m^  is  unable  to  adequately  maintain  its  position,  ?,  from  the 
reference  orbit  at  an  altitude,  h,  for  a  mission  time,  life,  using  a  thrusters  with  specific 
impulse,  Isp  and  efficiency  t|. 


The  values  for  the  three  constraints  listed  above  can  be  found  as  follows: 


1.  Fora  given  cluster  size  and  altitude,  the  accelerations  along  y  can  be  found  by 
evaluating  equation  (3-4) 

2.  The  corresponding  thrust  to  counteract  those  accelerations  is  given  by 

Ftot=ms/cr  (3-18) 

3.  Knowing  the  acceleration  and  assuming  continuous  thrusting,  the  total  AV  required 
is 


AVtot  =  drx|+|ry|+|rz|)(iife) 

4.  The  total  propellant  used  during  the  life  of  the  mission  is  simply 

mP4-e-^)ms/c 

5.  The  diameter  of  the  spacecraft  is 


(3-19) 


(3-20) 


Ds/c  ~ 


6  ms/c  ]3 


\n  Ps/c 

6.  Using  mp,  the  diameter  of  the  propellant  tank  is 


(3-21) 


D 


tank  — 


6  nip 

*  Pp 


7.  The  total  power  available  to  the  spacecraft  is 

PS/c  =  ms/cas/c 


(3-22) 


(3-23) 
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8.  For  electric  propulsion,  a  power  plant  is  required  on-board  the  satellite  to  provide 
power  to  the  thruster.  Therefore,  the  maximum  required  input  power  to  the 
thrusters,  using  the  current  value  of  Isp  being  evaluated  and  thrust  levels  from  step 
2  is 

Pb  =  jglfrgg)  (3.24) 

2r\ 

9.  The  mass  of  the  power  plant  can  be  found  using  input  power  calculated  in  step  8 
along  with  the  specific  power  of  the  power  plant 

mpp  (3-25) 

«PP 

Table  3-3  lists  the  maximum  values  chosen  for  these  constraints  as  well  as  other  input 
parameters. 


Table  3-3  Baseline  parameters 


“nn . 

lOW/kg 

^s/c 

2  W/kg 

Pd 

500  kg/m3 

Ps/c 

79  kg/m3  | 

life 

5  yrs 

ms/c 

100  kg 

(n^p/rris/c)max 

0.1 

(Pir/P  s/c)max 

0.20 

(^tank/^  s/c)max 

0.33 

1  [(Md+M  J/MJmax  0.30  ] 

3.3.3  Analysis  of  Cluster  Missions 

The  methodology  developed  in  the  previous  section  can  now  be  used  to  determine  the 
range  of  thruster  specific  impulse  and  efficiency  which  can  maintain  the  relative 
positions  of  a  satellite  in  a  cluster  without  exceeding  the  satellite  design  ratios.  Figure 
3-11  through  Figure  3-13  illustrate  the  feasible  regions  of  thruster  specific  impulse  vs. 
efficiency  for  clusters  orbiting  at  various  altitudes.  At  specific  impulses  below  the 
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feasible  regions,  the  required  propellant  to  maintain  the  cluster  during  the  five  year 
mission  life  exceeds  the  assumed  10%  maximum  initial  propellant  mass  fraction.  In 
actuality,  the  thruster  must  operate  at  an  1^  above  this  lower  limit  since  propellant 
utilization,  t|u=1,  is  assumed  in  this  initial  analysis.  Combinations  of  specific  impulse 

and  efficiency  above  and  to  the  left  of  the  feasible  regions  result  in  thruster  power 
requirements  greater  than  the  allowable  20%  of  spacecraft  power.  Because  the 
thruster  power  is  determined  by  equation  (3-24)  an  increase  in  thruster  efficiency 
allows  for  higher  specific  impulses  without  exceeding  the  power  limitations. 

As  shown  in  Figure  3-1 1  a  satellite  in  a  cluster  orbiting  at  1000  km  altitude  with  a  25  m 
baseline  requires  thrusters  operating  at  a  minimum  specific  impulse  and  efficiency  of 
2000s  and  40%  respectively.  SPT’s  and  ion  engine  technology  can  currently  achieve 
these  ranges.  If  the  baseline  is  increase  to  35  m,  thruster  requirements  increase  to 
2700s  specific  impulse  and  80%  efficiency.  Similar  effects  can  be  seen  for  clusters 
orbiting  at  10,000  km  (Figure  3-12)  and  at  GEO  (Figure  3-13).  At  10,000  km,  cluster 
baselines  in  the  hundreds  of  meters  can  be  achieved  using  SPT  or  ion  engine 
technologies.  These  baselines  increase  to  4000-6000  m  for  clusters  stationed  at  GEO. 
At  no  altitude  can  chemical  propulsion  maintain  cluster  formations  for  the  5  year  life. 
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Figure  3-1 1  Feasible  Isp  vs.  r|  at  1 ,000  km 
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Efficiency 


Figure  3-12  Feasible  Isp  vs.  r\  at  10,000  km 


Efficiency 

Figure  3-13  Feasible  Isp  vs.  rj  at  GEO 

Moving  a  cluster  to  higher  orbits  dramatically  increases  the  allowable  baselines  for 
sparse  apertures.  However,  due  to  the  increased  range  to  target  from  GEO,  ground 
resolutions  increase  by  only  about  a  factor  of  5  (see  Table  3-4).  The  achievable 
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resolutions  from  the  sparse  aperture  baselines  shown  in  Figure  3-11  through  Figure  3- 
13  are  listed  in  Table  3-4  and  Table  3-5. 


The  achievable  resolutions  from  the  sparse  aperture  baselines  shown  in  Figure  3-11 
through  Figure  3-13  are  listed  in  Table  3-4  and  Table  3-5. 


Table  3-4  Feasible  Sparse  Aperture  Ground  Resolution 


ALTITUDE 

CLUSTER 

BASELINE 

m 

HI 

HI 

1000  km 

25  m 

20  mm 

40  cm 

2km  | 

35  m 

14  mm 

■KK3I 

1 .5  km  I 

10,000  km 

25  m 

200  mm 

400  cm 

20  km 

100  m 

50  mm 

100  cm 

5km 

200  m 

25  mm 

50  cm 

2.5  km 

400  m 

13  mm 

25  cm 

1.3  km 

35768  km 

25  m 

710  mm 

15  m 

72  km 

2000  m 

9mm 

18  cm 

900  m 

5000  m 

4  mm 

7cm 

350  m 

7000  m 

3  mm 

5cm 

250  m 

For  mission  where  angular  resolution  is  the  main  concern  (such  as  for  imaging  of 
astronomical  objects),  Table  3-5  shows  that  angular  resolution  can  be  increased  by  a 
factor  of  order  1000  if  the  cluster  is  positioned  in  GEO  rather  than  LEO.  These  values 
can  be  compared  with  the  angular  resolutions  given  for  terrestrial  interferometers  in 
Table  3-2. 


Table  3-5  Achievable  Sparse  Aperture  Angular  Resolutions 


ALTITUDE 

CLUSTER 

BASELINE 

VISIBLE 

(.5mm) 

IR 

(10mm) 

RF 

(5cm) 

1000  km 

25  m 

4x1  O'3 " 

82x1  O'3" 

412" 

35m 

3x1 03" 

59x1 0'3" 

294" 

10,000  km 

25  m 

4x1  O’3" 

82x1  O'3" 

% 

CVI 

5 

100  m 

1x1  O'3" 

20x1  O'3" 

103" 

250  m 

41 2x1  O'6" 

8x10 3" 

40" 

400  m 

258x10-®" 

5x1  O'3" 

25" 

35768  km 

25  m 

4x1  O'3" 

82x1  O'3" 

412" 

2000  m 

51x1 0'6" 

Ixio-3" 

5" 

5000  m 

20x1  O'6" 

415x10-®" 

2" 

7000  m 

15x1 0'6" 

294x10-®" 

1.5" 
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Figure  3-14  shows  the  thrust  levels  per  unit  of  spacecraft  mass  required  by  a  satellite 
displaced  from  the  reference  orbit  by  15m,  150m,  and  2750m  (30  m,  300m,  5500m 
cluster  baselines  respectively).  These  values  were  found  using  the  full  length  equation 
(3-5).  Combinations  of  cluster  baselines  and  altitudes  which  can  be  achieved  with 
moderate  propulsive  requirements  (Figure  3-11  through  Figure  3-13)  are  all  seen  from 
Figure  15  to  require  thrust  levels  on  the  order  of  15pN  per  kg  of  spacecraft  mass  along 
y.  Thrust/mass  along  z  is  approximately  5  orders  of  magnitude  smaller,  justifying  the 
use  of  the  linearized  equation  (3-4).  In  addition,  thrust  levels  along  x  equal  zero. 

This  is  also  apparent  in  Figure  3-15  which  shows  the  required  AV  for  a  five  year 
mission.  The  AV  expended  in  z  is  negligible  compared  to  that  expended  in  y. 


Figure  3-14  Thrust/mass  Levels 
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Figure  3-15  Required  AV 


3.4  SUMMARY  AND  CONCLUSIONS 

This  chapter  has  examined  the  propulsion  system  requirements  for  maintaining  a  local 
satellite  cluster  formation  in  Earth’s  orbit  for  the  purpose  of  forming  sparse  aperture 
arrays.  Near  continuous  thrusting  by  the  propulsive  system  was  found  necessary  to 
maintain  the  relative  positions  of  the  satellites  within  allowable  tolerances.  Satellite 
mass,  volume,  and  power  constraints  limit  reasonable  cluster  baselines  to 
approximately  30  m,  300  m,  and  5000  m  at  1000  km,  10,000  km,  and  GEO  altitudes 
respectively.  To  maintain  these  cluster  baselines,  the  propulsive  system  must  operate 
at  minimum  1^  and  efficiency  equal  to  approximately  3000s  and  65%  respectively  with 
a  thrust/spacecraft  mass  ratio  of  approximately  15pN/kg. 
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Chapter  4 

Modeling  Ion  Engine  Performance 


In  Chapter  3,  the  propulsion  requirements  for  maintaining  the  relative  positions 
of  satellites  in  a  local  cluster  were  determined.  To  maintain  reasonably  sized 
formations,  it  was  found  that  station-keeping  thrusters  required  specific  impulse, 
thruster  efficiency  and  thrust  level  characteristics  similar  to  those  of  ion 
thrusters.  The  rest  of  this  thesis  examines  the  performance  of  micro  ion  engines. 
This  chapter  formulates  an  analytical  model,  known  as  Brophy’s  Theory,  to 
predict  the  performance  of  cylindrical  ring-cusped  ion  engines.  In  later 
chapters,  Brophy’s  model  will  be  used  to  predict  the  performance  of  cylindrical 
ion  engines  as  they  are  scaled  down  in  size  three  orders  of  magnitude  and  a 
linear  ion  microthruster  concept. 


4.1  OPERATION  OF  ION  ENGINES 

Figure  4-1  shows  the  cross  section  of  a  typical  cylindrical  ion  engine.  A  neutral 
gas,  such  as  Xenon  or  Argon,  is  fed  into  the  discharge  chamber.  Electrons  are 
injected  into  the  chamber  from  a  cathode  emitter.  On  occasion,  the  electrons 
collide  with  the  neutral  atoms  with  sufficient  energy  to  ionize  one  neutral 
resulting  in  a  positively  charged  ion  and  one  or  more  additional  electrons. 
During  steady  state  operation,  the  neutrals,  ions  and  electrons  form  a  quasi¬ 
neutral  plasma.  Ions  in  the  vicinity  of  the  grid  are  accelerated  across  a  potential 
difference  forming  an  ion  beam  exiting  the  chamber  and  producing  thrust.  In 
some  configurations  a  third  grid  surface  is  added  to  control  ion  velocity  (specific 
impulse)  without  affecting  ion  mass  flow  rate  (thrust).  Electrons  in  the  chamber 
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eventually  reach  the  anode  where  they  either  complete  the  circuit  back  to  the 
cathode  or  are  injected  into  the  ion  beam  to  neutralize  it. 

VNET 


4 


^acc 


VTOt 

Figure  4-1  Current  Balances  in  Ion  Engines 


The  following  is  a  list  of  definitions  for  the  terms  shown  in  Figure  4-1 

JB=  Beam  ion  (and  neutralizer  electron  current) 

JE  =  cathode  emitted  current 

Jc  =  ion  current  to  cathode-potential  surfaces 

JD  =  current  through  discharge  power  supply 

Jp  =  total  ion  production  rate 

Ja  =  Ion  current  to  anode 

J**  =  Ion  current  intercepted  by  accelerator  grid 

VD  =  Discharge  potential  (potential  difference  between  cathode  and  plasma) 

V^,  =  Neutralizer  potential  (potential  difference  between  plasma  and  ion  beam) 
VB,  =  Acceleration  grid  potential 

4.2  BROPHY’S  MODEL 

Brophy’s  model  [1]  uses  the  conservation  equations  for  mass,  charge,  and 
energy  to  predict  the  beam  ion  production  cost,  eB  is  the  energy  necessary  to 

produce  an  ion  which  becomes  part  of  the  beam.  The  derivations  in  this 
chapter  are  a  more  annotated  version  of  those  found  in  reference  [2]. 
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4.2.1  Fundamental  Equations 

Beam  ion  production  energy  can  be  defined  as  the  total  energy  expended  per 
beam  ion  less  the  useful  energy  per  beam  ion.  Equations  (4-1)  and  (4-2)  define 
useful  power  and  total  power  respectively 


UsefiilPower  =  Jb(VB  +  Vp  )  (4-1 ) 

Total  Power  =  J BVB  +  JD  VD  +  J ^  VB  +  Pheaters  (4-2) 

Using  equations  (4-1)  and  (4-2),  eB  can  be  found  as  shown  in  equation  (4-3). 

_  TotalPower - UsefulPower  _  (Jp  - Jb)^D  +  Jacc^B  +  Ph  /a 

-  T  —  Y  '  ' 

JB  JB 

Examining  Figure  4-1  the  following  current  balance  can  be  found 

Jp-JB  =  JE  +  Jc  +  Jacc  (4-4) 


Substituting  equation  (4-4)  into  equation(4-3)  and  expanding  yields 

eB  = 


+  jp  +  JCVD  +  Jacc(VB  +Vd) 
V  JP  JJB  JB 


This  can  be  rewritten  as 


eB=^  +  ^VD+^(VB  +  VD) 
IB  IB  iB 


where 


£p  - 


JEVD  +PH 


(4-5) 


(4-6) 


(4-7) 


and  %  =  ■£.%=£.  4*  =7* 

Jp  Jp  Jp 
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4.2.2  Discharge  Energy  Balance 

Equation  (4-6)  defines  in  general  terms  the  expression  for  beam  ion  production 
cost  The  goal  from  here  on  is  to  express  the  terms  in  equation  (4-6)  as 
parameters  which  are  known  or  can  be  measured.  First,  looking  at  the  terms  in 
the  equation  for  Ep  an  energy  balance  for  the  discharge  chamber  can  be 
expressed  as 


u++EJj 


U: 


j  |  JLPVD 


h£B 


Op  4- Jr  ~  Jja  ~JlP^ 


JPV, 


EVD  _ 


£P 


fH 

Jp 


(4-8) 


where 

U+  =  ionization  energy  per  ion 
Uj  =  excitation  energy  of  level  j 
Jj  =  excitation  rate  (to  level  j) 

Jlp  =  loss  rate  of  primary  electrons 

em  =  mean  energy  of  Maxwellian  electron  group 


Because  the  first  two  terms  in  equation  (4-8)  are  known  constants  specific  to  the 
type  of  propellant  we  can  define 


U5 


(4-9) 


Substituting  equation  (4-9)  into  equation  (4-8)  and  expanding  gives 

JLPJE  Jia 


^^  =  e0  +  em+(VD-em) 

J  n 


JEJp  Jp 


em  em 


f  T  > 
JE 


vJp; 


Solving  for  yields 


JEVD 


£0+£ 


m 


f  J.  ^ 
1-lia 


V 


) 


1- 


Jlp 

JE 


\(  9  \ 

1 


A 


VD 


(4-10) 


(4-11) 


Substituting  equation  (4-11)  into  (4-7)  and  writing  PH  =  JEVC,  the  equation  for  £p 
becomes 
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(4-12) 


We  can  define 


eP 


"f"  em 


h  ) 


»LP 

je  ; 


V 


1  em 


'Dj 


1  + 


V| 

vdJ 


[£0+Em(f  fja)] 


1  +  - 


*Dj 


so  that  equation  (4-12)  becomes 


1- 


JLP 
JE  ) 


(4-13) 


(4-14) 


£p*  is  seen  to  be  the  energy  per  ion  created  if  no  primary  electrons  were  to 
escape.  The  terms  for  ep*  in  equation  (4-13)  will  be  related  later  in  equation  (4- 
32)  to  known  or  measurable  parameters.  Now,  however,  let  us  express  JLP/JE 
in  known  or  measurable  parameters. 


4.2.3  Survival  Equation  for  Primary  Electrons 

Jlp/je  =e_CT‘°*n”1'  (4-15) 

where  atot=  (cr+  +aexc)primaijes. 


Also,  le  is  the  path  for  a  primary  electron  before  it  would  be  captured  by  the 
anode,  if  it  did  not  collide  with  a  neutral  before  that.  This  path  length  is  that  of 
the  electron’s  helical  path  around  one  of  the  magnetic  lines  of  force  created  by 
the  confinement  magnets. 

The  neutral  density  in  equation  (4-15)  is  related  to  the  flow  rate  by 
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The  quantity  C0  is  a  measure  of  the  confinement  effectiveness  for  primary 
electrons  (better  for  long  electron  path  le,  small  grid  open  area  If  C0->~, 

the  energy  cost  per  ion,  £p,  tends  to  the  limit  e^,  which  then  represents  the  cost 
per  ion  with  no  primary  losses. 


Keep  in  mind  that  the  ultimate  goal  is  to  find  equation  (4-6),  in  terms  of 
known  or  measurable  parameters.  Up  to  this  point,  we  have  been  working  to 
express  ep  in  such  terms.  Equation  (4-18)  represents  the  progress  thus  far,  but 
it  still  remains  to  express  e,,  in  such  terms.  This  is  the  focus  of  the  next  section. 
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4.2.4  Calculation  of  Primary  /  Secondary  Population  Ratio 

Primary  electrons  are  endowed  initially  with  an  energy  VD,  and,  if  they  did  not 
escape,  would  all  thermalize  eventually,  to  an  energy  ^  The  energy 
transferred  to  the  neutrals  by  each  primary  is  (VD  -  ej.  In  steady  state  operation, 
this  energy  loss  by  the  primaries  must  exactly  balance  the  energy  gained  by  the 
neutral  propellant.  Transfer  of  energy  from  the  primary  electrons  is  assumed  to 
occur  in  four  ways;  1)  ionization  of  atoms  by  primaries,  2)  excitation  of  atoms  by 
primaries,  3)  ionization  of  atoms  by  secondaries  (Maxwellian),  4)  excitation  of 
atoms  by  secondaries.  Thus  the  balance  of  energy  rate  per  unit  volume  is  , 

^int  o  =  ^pi  "*■  ^pe  ^si  ^se  (4-1 9) 

The  rate  at  which  primaries  disappear  (i.e.  thermalize  to  secondaries)  is  simply 
the  rate  of  ionization  or  excitation  by  primaries  (a  primary  is  assumed  to  become 
a  secondary-Maxwellian  after  one  ionization  or  one  excitation).  So,  the  net 
energy  input  rate  per  unit  volume  due  to  injection  of  primaries  is  (assuming  the 
primaries  do  not  leave  the  chamber  without  ionizing  or  exciting  a  neutral) 

^into  =  nnI1pvp<ytot(^D)(^D  — em)  (4-20) 

where  ot0t=  (o+  +cexe) 

As  indicated  in  equation  (4-19),  this  energy  is  used  by  the  primaries  and  their 
secondary  “progeny”  to 

(1)  Produce  ionization  by  the  primaries.  Per  ionization  event,  this  uses  U+ 
+  em,  since  the  new  electron  created  has  energy  em.  The  net  energy  input 
rate  per  unit  volume  due  to  ionization  by  primaries 

Ep,+  =  nnnpvpo+(VD)(u+  +  em)  (4-21) 

(2)  Excite  atoms  by  primaries.  Total  energy  rate  per  unit  volume  is 

^p,exc  =  nnnpvp°exc(^D)Uexc  (4-22) 
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(a  shorthand  for  nnn  pVpXojUj) 

j 

(3)  Produce  ionization  by  secondaries  (Maxwellian).  The  rate  per  unit 
volume  is 

oo 

Esi  =nn  Jfm  (c)co+  (c)4jic2dc  a  nncenmc+  (4-23) 

0 

1  00 

where  a+  =  — - f  fm  (c)ca+  (c)4rcc2dc 

cenm  J 

oo 

.  nm  =  Jfm47tc2dc 
0 

The  ionization  cross-  section  o+(c)  is  zero  below  c+  =  ^2eU+/me .  Using  a 
Maxwellian  form  for  fm(c),  we  find 

oo 

otot  =  Je“uua+(u)dii  (4-24) 

u* 

/  E 

where  u  = - 

l  kTe 

and  the  energy  spent  by  secondaries  in  ionization  (per  unit  time  and  volume)  is 
then 

Em,+  =nnnmcec+(u++em) 

Similarly,  the  energy  spent  in  excitation  is 

Em,exc  =  nnnmce<5excUexc 


(4-25) 

(4-26) 
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npvp^tot(^D  ^m)  npvp^+(^D)(^  +  ^ra)  +  (7exc(^D  )^exc  j  +  nmce  [^+(^  +^m)  +  ^cxc^exc  j 


This  can  be  solved  for 


lp  . 


‘m 


n 


lT+£m+U 


'exc 


P  _  ce 


m  1  wexc  — 


n®  VP  (VD  -  £m )  gtol( Vp )  -  fu+  +£m)-gexi(Vp) 

0+  C+  v  '  G+ 


(4-27) 


Equation  (4-27)  is  important  because  it  expresses  the  ratio  of  primary  to 
secondary  electrons  in  the  chamber  in  terms  of  known  or  measurable 
parameters.  Note  that  it  is  a  function  of  Te  for  a  fixed  VD. 


Before  continuing  on  with  the  derivation  of  £„',  we  must  return  to  the  expression 
for  e0,  equation  (4-9).  Recalling  the  shorthand  used  in  equation  (4-22)  to 
represent  the  sum  of  energy  levels,  equation  (4-9)  can  be  written  as 

J, 


£0  =U+  + 


'exc 


U 


exc 


(4-28) 


Utilizing  the  rate  per  unit  time  and  volume  terms  from  equations  (4-21),  (4-22)  , 
(4-25),  and  (4-26),  the  quantity  Jexc/Jp  in  equation  (4-28)  is 


JeXc  _  Jp.exc  ^m.exc  _  npvp°exc(^D)^' nmce<Jexc 


Jp  Jp,+  +Jm,+  npvp<y+(^D)  + nmcecy+ 


(4-29) 


Do  not  confuse  the  terminology  used  here.  Jp  represents  the  current  of  ions 
produced  (through  ionization  by  either  by  primary  or  secondary  electrons)  while 
Jp  +  represents  the  part  of  Jp  produced  by  primary  electrons. 


Substituting  equation  (4-29)  into  equation  (4-28)  yields 


_TT+  TT  npvpcrexc(^D)  +  nmce<yexc 

£0=  6XC  nP  vptf  +  (VD )  +  nmcea+ 


(4-30) 
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and  this  does  depend  on  Te  and  — (Te).  Substituting  equation  (4-27)  for 

nm 

puts  e0  in  terms  of  known  or  measurable  parameters. 


Having  found  an  expression  for  e0,  we  can  continue  with  determining  £p\  The 
expression  for£p“,  equation  (4-13),  can  be  simplified  by  neglecting  ion  capture 
by  the  screen,  fja=0,  and  heating  power,  Vc=0  so  that 


*  _  e0  ~*~^m 


£P  “ 


j_£m_ 

VD 


(4-31) 


The  expression  fore0,  equation  (4-30),  can  be  substituted  into  equation  (4-31). 

°D 

utilizing  the  expression  for  — -,  equation  (4-27),  and  simplifying,  we  obtain 


ep*  =  VDa+(VD) 


a+(u++£m)  +  CTexcUexc 


[^exc^+  (VD ) - Cexc ( VD )CT+  ]u exc  +  a+  ctot (VD )(VD  “Em) 


(4-32) 


NOTE:  An  intermediate  expression  for  Ep*  (still  containing  — — ),  which  will  be 


*m 


useful  later,  is 


'D 


* 

eP  = 


'  tot 


V  a+ 


j  I  nm  ce 


(4-33) 


nP  VP 


^+(vd) 


At  this  point  we  have  found  an  expression  for  £p*  in  terms  of  measurable 
parameters  for  use  in  defining  the  expression  for  e^,,  equation  (4-18).  The 
expression  for  £p  is  to  be  used  for  defining  £g,  equation  (4-6).  What  remains  to 
be  found,  however,  is  in  expression  for  propellant  utilization,  t|u,  present  in 
equation  (4-18). 
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4.2.5  Calculation  of  Utilization  Efficiency 

The  derivation  of  the  expression  for  propellant  utilization,  tju,  begins  with  a 

balance  of  the  rate  ions  are  produced.  Specifically,  the  number  of  ions 
produced  per  second  equals  the  rate  at  which  primaries  create  ions  plus  the 
rate  at  which  secondaries  create  ions.  Using  the  ionization  rates  per  unit 
volume  from  equations  (4-21)  and  (4-25),  this  balance  becomes 


IP 

e 


=  (vol)(np' vpa+  (VD)  +  nmcec+  )nn 


where  nn  =  m(l-Tlu)— — - 

c^nAgmi 

Noticing  that 


(4-34) 


Jr>  =  fpJn  = 


Rj  p  “ 

r  m  . 


m; 


(4-35) 


expands  equation  (4-34)  to 

(vol)(np  vpa+  (VD )  +  nmceo+ =  7^ 

cnYnAg  rB 

But  also  n+  =  nm  +  np  and  JB  =  (0.61)en+vBAg<(>i  so 


(4-36) 


n 


m 


p  (0.61)evBAg<|>j 


(4-37) 


Dividing  the  left  hand  side  of  equation  (4-36)  by  the  left  hand  side  of  equation 
(4-37)  and  the  right  hand  side  of  equation  (4-36)  by  the  right  hand  side  of 
equation  (4-37)  and  simplifying  yields 


1  + 


np  vp(Vd)c+(Vd) 


(vol)ceo+- 


lm 


G+  4(1  -T1U)  (O^OevEAg^i 


i  +  - 


n. 


Cn^nAg 


JBfB 


nu 


n 


m 


Equation  (4-38)  can  be  rearranged  so  that 


(4-38) 
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l-'Hu 

nu 


(0.6  l)ev  B  Ag24>n<t>|Cn 


1  +  - 


*m  J 


4(vol)JBfBceC+ 


1  + 


np  vp(vd )q+(vD)>| 

J 


(4-39) 


nm  ce 


Notice  that  the  bracketed  term  in  the  denominator  of  equation  (4-39)  appears  in 
equation  (4-33).  Swapping  out  that  term  results  in 


l-'Hu  _ 

nu 


(0.61)evBAg2<t>n<t>iCn 


1  + 


np' 


n 


m  ) 


4(vol)JBfBceo+ 


^VD  atot(VD)  nP  VP^ 


vep* 


nm  ce 


(4-40) 


We  can  solve  equation  (4-40)  for  r|u  so  that 

1 


Tlu  = 


1+y 


(4-41) 


where  y  = 


0.61 


8p*evBAg2$n$iCn 


1  + 


n. 


) 


(vol)JBfBatot(VD) 


n. 


Vnm 


Figure  4-2  summarizes  the  procedure  for  using  Brophy’s  model. 
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4.3  ARAKAWA  ALGORITHM 


In  the  previous  section,  Brophy’s  model  was  derived  allowing  the  energy  costs 
per  beam  ion,  to  be  calculated  from  a  simple  algebraic  equation.  Although 
gas  properties,  chamber  geometry,  and  operating  conditions  are  known  or 
easily  found,  several  important  parameters  must  still  be  determined. 
Specifically,  the  pathlength  of  an  electron,  le,  is  dependent  on  the  magnetic  field 
within  the  chamber.  Also,  the  fraction  of  ions  reaching  the  beam  and  cathode 
potential  surfaces,  fB  and  fc  respectively,  are  dependent  on  the  plasma  density 
at  the  surfaces.  Plasma  density,  in  turn,  is  partially  dependent  on  diffusion  of 
electrons  across  and  along  magnetic  field  lines.  Arakawa’s  algorithm  [3]  was 
developed  to  calculate  the  magnetic  field  in  the  chamber  and  implement 
Brophy’s  model. 

4.3.1  Magnetic  Field  Analysis 

From  Ampere’s  law  we  know  that 

VxB=p0J  (4-43) 

which  for  the  case  of  no  current  source,  J,  this  can  be  expanded  to 

VxB=Vx(VxA)  =  V(V  •  A)  -  V2  A  =  0  (4-44) 

where  A  is  of  the  magnetic  potential  vector.  Given  that  the  magnetic  potential  is 
divergenceless  (i.e.  V  •  A  =  0),  this  reduces  to 


V2A=0  (4-45) 

The  location  and  direction  (N-S)  of  each  magnet  in  the  chamber  is  specified  in 
an  input  file  to  the  magnetic  analysis  code.  The  code  divides  the  chamber  into 
numerous  triangular  elements.  The  finite  element  method  is  applied  to  this  grid 
in  order  to  solve  for  the  magnetic  potential  at  each  node.  The  magnetic  field 
vector  in  cylindrical  coordinates  can  be  written  as 


B  =  Vx A  = 


ld\ 
r  ae 


dAe  V  f  r  3Ar  3AZ 

dz  J  v  3z  3r 
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which,  for  the  axisymetric  case,  reduces  to 


B  =  -^f+ 
dz 


r[^(rAe)]S  =  "r[^lA^]f  +  S(rA«))  <4'46) 


For  convenience  then,  the  magnetic  analysis  code  outputs  to  a  file  the  value  of 
rAe  at  each  node  location. 

4.3.2  Primary  Electron  Containment  Length 

Having  calculated  the  value  of  rAg  at  each  node,  a  Monte  Carlo  simulation  of  the 


electrons  in  the  chamber  is  carried  out.  Starting  from  the  cathode  location,  an 
electron  is  given  an  energy  corresponding  to  the  discharge  voltage.  Its  position 
and  velocity  are  surveyed  as  a  function  of  time  during  its  path  through  the 
magnetic  field.  Its  trajectory  is  computed  using  the  Runge  Kutta  method 
according  to  the  following  equations  of  motion  while  preserving  angular 
momentum,  M0. 


dvz  _  rv  ld(rA0) 
dt  0  r  3r 

(4-47) 

dV'=meVe2  eveiS(rAe) 

dt  e  r  r  dr 

(4-48) 

=  merv0  -ctAq  =  constant 

(4-49) 

While  in  the  chamber,  a  probability  exists  that  the  electron  will  have  its  trajectory 
altered  either  by  an  elastic  collision  with  a  neutral  atom  or  by  the  “wave”  of  the 
oscillating  plasma.  These  probabilities  are  calculated  from  the  mean  free  paths 
of  the  respective  phenomena.  The  electron  continues  until  it  reaches  the  anode 
location.  The  path  lengths  of  the  several  thousand  electrons  simulated  are 
averaged  to  yield  the  primary  electron  confinement  length,  le  needed  in 
Brophy’s  model.  The  time  the  electrons  spend  in  each  triangular  grid  element  is 
totaled  to  provide  the  distribution  of  primary  electrons  within  the  chamber. 
Because  the  simulation  is  estimating  the  electron  path  length  assuming 
ionization  does  not  occur,  no  secondary  electrons  are  created  (see  section 
4.2.3). 
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4.3.3  Discharge  Chamber  Performance 

Having  found  the  electron  densities,  the  Arakawa  algorithm  proceeds  to 
determine  the  plasma  distribution  in  the  chamber.  The  governing  equation  for 
the  particle  balance  within  the  chamber  is 


-V-([D][Vn])  =  Q  (4-50) 

where  Ddescribes  the  spatial  variation  of  the  diffusion  coefficients  parallel  and 
normal  to  the  magnetic  field  lines  and  Q  is  the  ion  production  rate  per  unit 
volume.  Reference  [3]  provides  a  fuller  development  of  those  terms.  The  result 
is  the  differential  equation 


1  a/'  9n 

r3r(.  2,1 


3n 

dz 


dn 


dn) 


+£l>§ +D"£r6=o  <4-5i) 


where 


Dj  j  =  Dp  cos2  9  +  sin2  9 

D22  =  Dp  sin2  9  +  Dn  cos2  9 
D12  =  D21  =  (Dp  —  Dn)sui9cos9 

and  Dp  and  DN  are  diffusion  coefficients  parallel  and  normal  to  magnetic  field 
lines,  respectively.  Plasma  density  within  the  chamber  is  found  by  solving 
equation  (4-51). 


The  ion  flux  to  surfaces  is  assumed  to  be  proportional  to  plasma  density  at  that 
area.  Knowing  the  plasma  density  along  the  exit  and  cathode  potential 
surfaces,  the  fraction  of  ions  to  these  regions,  fB  and  fc  respectively,  can  be 
calculated.  The  values  for  le,  fB,  fc  are  used  together  with  other  chosen  or 
constant  values  in  Brophy’s  model  to  calculate  the  energy  cost  per  beam  ion,  eB. 


Figure  4-3  summarizes  Arakawa’s  algorithm. 
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Input  Data 


Figure  4-3  Summary  of  Arakawa  Algorithm 

4.4  SUMMARY 

In  this  chapter,  the  model  used  to  predict  the  performance  of  ring-cusped  ion 
engines  was  derived.  It  was  shown  that  energy  cost  per  beam  ion,  %,  can  be 

formulated  as  a  simple  algebraic  equation  known  as  Brophy’s  model.  Most  of 
the  parameters  in  Brophy’s  model  are  known  a  priori  given  the  choice  of 
propellant  gas,  chamber  geometry,  and  operating  conditions.  The  remaining 
parameters  can  be  determined  using  Arakawa’s  algorithm  the  fundamentals  of 
which  were  also  described  in  this  chapter.  The  next  chapter  uses  this  model  to 
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examine  how  performance  of  cylindrical  ion  engines  changes  as  their  scale  is 
reduced  several  orders  of  magnitude. 
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Chapter  5 

Scaling  of  Cylindrical  Ion  Engines 


The  previous  chapter  formulated  Brophy’s  model  and  outlined  Arakawa’s  code  which 
implements  the  model.  The  code  is  discussed  first  with  an  emphasis  on  the 
formulations  used  to  simulate  particle  behavior  within  the  chamber.  This  is  important 
because  of  the  assumption  in  subsequent  analysis  that  the  physics  remain  the  same  at 
smaller  scales.  Next,  an  argument  is  outlined  as  to  how  thruster  parameters  such  as 
mass  flow  rate  and  discharge  voltage  should  change  as  thruster  size  is  scaled.  The 
effect  on  performance  of  scaled  cylindrical  ion  thrusters  is  then  presented. 


5.1  SCALING  ARGUMENT 

This  section  examines  how  each  of  the  input  parameters  to  the  Arakawa  code  should 
scale  in  order  to  maintain  the  same  discharge  chamber  performance  as  chamber  size 
is  reduced.  As  developed  in  Chapter  4,  the  average  beam  ion  energy  cost  is  given  as 

£b=-t - A— - nr+fvD  (S'1) 

fB{l-exp[-C0m(l-iiu)J}  fB 

where  C0  =  4a0le/ev0Agtp0  and  m  is  the  mass  flow  rate  per  chamber  in  Aeq/s.  Equation 

(5-1)  indicates  that  it  is  desirable  to  maximize  the  value  of  the  primary  electron 
containment  length,  le,  in  order  to  reduce  the  power  needed  to  create  beam  ions.  In 
other  words,  the  longer  an  average  primary  electron  remains  in  the  chamber  before 
being  lost  to  the  anode,  the  more  effective  ionizer  it  is  and  the  more  efficient  the 
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thruster  will  be.  Ring  and  line  cusped  ion  engines  increase  le  by  employing  magnets 
which  cause  the  electrons  to  spiral  around  magnetic  field  lines. 

Equation  (5-1)  also  indicates  that  it  is  desirable  to  maximize  the  fraction  of  ions 
reaching  the  beam,  fB,  and  minimize  the  fraction  of  ions  reaching  cathode  potential 
surfaces,  fc.  In  a  well-designed  chamber,  the  presence  of  a  magnetic  field  tends  to 
contain  the  primary  electrons  near  the  center  of  the  chamber.  Because  a  higher 
percentage  of  ionizations  occur  near  the  chamber  center  and  the  ion’s  ability  to  diffuse 
across  field  lines  is  hindered,  fB  increases  and  fc  decreases  yielding  improved 
performance. 

As  the  size  of  the  thruster  is  decreased,  the  chamber  exit  area  will  scale  as  Ag~L2 
where  L  is  a  characteristic  length.  To  maintain  the  same  energy  cost  per  beam  ion  as 
chamber  size  is  scaled,  equation  (5-2)  suggests  that  the  primary  electron  confinement 
length  and  mass  flow  rate  should  scale  as  le~L  and  m~L. 


Maintaining  le  proportional  to  a  characteristic  chamber  length  should  be  a  matter  of 
keeping  the  same  “degree  of  particle  interaction”  between  the  electrons,  neutrals,  and 
ions.  Ideally  this  would  be  accomplished  by  altering  the  magnetic  field  strength  and 
particle  densities  so  that  the  Larmor  radius  and  mean  free  paths  scale  as  L.  Magnetic 
field  strength,  B,  and  Larmor  radius  are  related  by 


(5-2) 


mevn 
B  =  — 
erL 

Equation  (5-2)  implies  that  B-1/L  since  me  and  e  are  constants  and  vp  is  proportional 
to VVd  which  is  assumed  invariant.  The  mean  free  paths  for  elastic  collisions  of 
electrons  with  neutrals  and  ions,  inelastic  ionizing  collisions  with  neutrals,  and 
electron-wave  interactions  due  to  plasma  oscillations  are,  respectively 


^e-n  — 


no^e-n 


(5-3) 
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(5-4) 
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^anom  ”  2tt| 


^£q  Vp  A 
V  e  n  +/ 


X 


=  271^ 


(5-6) 


where  n+=np+nm  since  a  primary  electron  is  assumed  to  become  a  secondary  electron 
after  colliding  with  a  neutral  to  form  an  ion.  For  the  mean  free  paths  to  scale  as  L, 
equations  (5-3)  through  (5-5)  require  that  the  neutral  and  ion  densities  scale  as  n  ~1/L 


and  n+~1/L.  With  n+~1/L  equation  (5-6)  indicates  that  A.anom  does  not  scale  as  L  but  as 

Vl  with  the  result  being  that  electron-wave  interactions  occur  less  frequently  as  size 
is  reduced.  Although  equation  (5-6)  represents  the  anomalous  mean  free  path 
used  in  Arakawa’s  primary  electron  trajectory  simulation,  it  may  not  be  appropriate 
under  the  assumption  under  Bohm  diffusion.  In  such  a  case,  Xanom  ~  r^^  =  rT^c/eB  so 


that  ~  L,  not  Vl  .  In  this  analysis  however,  equation  (5-6)  is  used. 


We  see  that  n0~1/L  is  consistent  with  m~  L  since 

_  4m(l--nu) 

mivnAg<|>o 

That  n+~1/L  is  consistent  with  m~L  since 


(5-7) 


JB  =mqu  =0.6l(n+vBAg<j>i)  (5-8) 

where  Ag~L2  and  vB  is  proportional  toJT^  which  is  assumed  invariant.  It  is  also 
desired  that  the  propellant  utilization,  T|u,  remain  invariant.  Propellant  utilization  can  be 
represented  as  formulated  in  Chapter  4 


nu  =1  / 


0.15e£p  vBAg  <|>i<t>ocn 


1  +  - 


1  +  - 


J 


JBfBVDvp(Vol)aT 


f  n  ^ 
_£ 

Vnm  ) 


(5-9) 
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where  V-^/Vd  ,ca~Jf^andi  and  so  remain  invariant  given  the  assumptions 

that  Te  and  VD  are  invariant.  Additionally,  Ag~L2,  JB~L  and  Vol~L3.  r|u,  then,  will  remain 
invariant  only  if  the  ratio  of  primary  to  secondary  electrons,  np/nm,  and  baseline  energy 
cost  per  ion,  ep*,  are  invariant.  np/nmand  ep”  can  be  calculated,  respectively,  as 

.  |  U++8m+Uexcf^l 
£e_  _ l  G+  J _ 

VpJ  (VD-em)f^]-f^V  +  em)-2s« 

l  va+y  va+/  '  °+ 

and 


Bp 


(5-11) 


Therefore,  np/nm  and  ep  will  remain  invariant  if  ce  is  invariant.  ce  is  given  by 


indicating  that  it  is  consistent  for  t|u  and  Te  to  simultaneously  remain  invariant  as 

thruster  scale  is  reduced.  A  more  general  scaling  argument  for  electric  thrusters  is 
given  by  Khayms  in  reference  [1]. 


A  10  cm  diameter  ion  thruster  typically  utilizes  permenant  magnets  which  produce  a 
0.1  T  magnetic  field  near  the  magnet  face.  Scaling  the  thruster  to  a  diameter  of 
approximately  1  cm  requires  a  1  T  magnetic  field  which  may  be  achieved  using  rare- 
earth  ceramic  magnets.  Achieving  the  10  T  field  strength  necessary  for  a  1  mm  size 
microthruster  would  presumably  require  superconducting  wire  solenoids  to  avoid 
ohmic  losses  (and  subsequent  thermal  issues)  in  the  wires. 
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5.2  SCALING  OF  CYLINDRICAL  CHAMBERS 


The  analytic  code  described  in  Chapter  4  was  used  to  determine  the  effect  that  scaling 
has  on  discharge  chamber  performance.  The  magnetic  analysis  part  of  the  code  was 
found  to  produce  erroneous  results.  In  accordance  with  Gauss’s  Law,  the  magnetic 
flux  through  a  closed  surface  should  equal  zero, 

<fB-dA  =  0  (5-13) 

However,  performing  this  check  with  the  field  generated  by  Arakawa’s  magnetic  field 
program  yielded  integrations  differing  significantly  from  zero.  This  indicates  a  flaw  in 
the  calculation  of  the  magnetic  potential.  Therefore,  a  commercial  package  known  as 
MAGNETO  [2]  was  used  in  its  place. 

A  7cm  ion  thruster  developed  at  Colorado  State  University  was  chosen  as  a  baseline 
because  it  has  undergone  experimental  testing  and  numerical  simulation  [3].  The 
chamber  geometry  with  corresponding  field  lines  and  flux  density  contours  calculated 
by  MAGNETO  are  shown  in  Figure  5-1  and 

Figure  5-2  respectively.  Magnetic  flux  density  was  held  at  1200  G  at  the  magnet 
surfaces  (representative  of  SmCo  magnets)  for  all  cases. 
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Figure  5-2.  7cm  CSU  flux  density  contours  (Gauss  units). 

The  structural,  material,  and  power  resources  required  to  produce  magnetic  field 
strengths  of  several  Tesla  in  the  chamber  of  an  ion  microthruster  are  not  yet  available. 
Therefore,  the  effect  of  scaling  a  chamber  while  holding  the  permanent  magnet 
surface  strength  constant  at  that  of  the  7cm  chamber  was  investigated.  Many  of  the 
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relationships  presented  in  the  previous  section  will  not  hold  since  the  magnetic  flux 
density  will  not  scale  proportional  to  1/L 

Figure  5-3  shows  beam  ion  production  cost,  sB,  as  a  function  of  propellant  utilization  for 

cylindrical  chamber  diameters  ranging  from  7cm  to  0.07cm.  Despite  best  efforts,  the 
values  for  the  7cm  case  differ  significantly  from  those  obtained  in  reference  [4].  The 
higher  eb  calculated  for  this  paper  are  due  to  a  lower  fraction  of  ions  reaching  the  beam 

(fB=0.21  vs.  0.29  in  Ref.  [3]).  This  discrepancy  is  probably  due  to  differences  in  the 
magnetic  field  as  generated  by  MAGNETO  and  Arakawa’s  magnetic  field  analysis 
program.  Arakawa’s  program  overstates  the  strength  of  the  magnetic  field  resulting  in 
a  higher  concentration  of  ions  near  the  chamber  center.  Consequently,  Arakawa’s 
program  predicts  a  larger  fraction  of  these  ions  reaching  the  beam.  Although  the 
values  of  eB  in  this  analysis  do  not  agree  with  experimental  data,  it  is  still  possible  to 

identify  the  factors  influencing  eB  as  chamber  size  is  scaled. 


0  02  0.4  0.6  0.8  1 

Propellant  Utilization 


Figure  5-3.  eBvs.  x\u. 

The  results  of  scaling  the  7cm  CSU  ion  thruster  as  generated  by  the  Arakawa  code  are 
listed  in  Table  5-1.  Mass  flow  rate  was  initially  set  at  0.15  Aeq  consistent  with  it’s 
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operation  at  CSU.  The  neutral  density  scales  as  1/L  as  expected  due  to  the  scaling  of 
mass  flow  rate  as  L.  Due  to  the  lack  of  a  magnetic  field  the  primary  electron 
confinement  length  scales  greater  than  L.  This  results  in  proportionately  less 
ionization  so  that  ion  density  scales  less  than  proportional  to  1/L.  The  equilibrium 
temperature  does  not,  therefore,  remain  invariant  but  increases  slightly  as  the  thruster 
is  scaled  down. 

Figure  5-4  shows  the  distribution  of  primary  electrons  in  the  7cm  chamber.  The  peak 
electron  density  occurs  at  the  position  of  the  cathode.  As  shown  in  the  figure,  the 
magnetic  field  tends  to  contain  the  primaiy  electrons  near  the  chamber  center. 
Ionizations  will  then  also  tend  to  occur  near  the  chamber  center  resulting  in  a  plasma 
density  distribution  as  shown  in  Figure  5-5.  A  large  fraction  of  ions,  fB,  reaches  the 
beam  tending  to  improve  performance.  Note  the  sharp  increase  in  eb  as  L  decreases 
indication  larger  wall  loses. 


Table  5-1.  Scaling  of  7cm  CSU  ion  thruster  (qu=0.7) 


units 

L 

0.1  L 

f  0.01L  1  Scales  as  1 

k _ 

m 

0.07 

0.007 

0.0007 

L 

KQSfSI 

0.15 

0.015 

0.0015 

L 

cm 

m 

1 .8e-4 

1.8e-4 

1.8e-4 

1 

EPH 

m 

1.01 

0.101 

0.0101 

L 

mm 

m 

334.0 

34.1 

3.44 

_L°" 

mm 

m 

4.7 

0.47 

0.047 

L 

mm 

mm 

9.0e-4 

2.8e-4 

09.1  e-5 

Vl 

m 

m 

2.09 

0.076 

0.00629 

wmssm 

m 

Hi! 

12.5 

46.4 

376.0 

k _ ! 

1018/m3 

4.76 

47.6 

476.0 

1/L 

EH 

lO^/m3 

0.130 

1.25 

12.7 

ssL"0'99 

EH: 

1018/m3 

0.0206 

0.310 

3.53 

=L'in 

Jb 

unitless 

0.21 

0.14 

0.12 

WEsami 

Jc  ....  . 

unitless 

0.79 

0.86 

0.88 

_|_-0.02 

i _ 

eV 

3.61 

3.78 

3.82 

srL”001 

W 

70.4 

72.5 

73.1 
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W 
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4029 
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Figure  5-4.  Primary  electron  densities  for  7cm  chamber. 


Figure  5-5.  Plasma  densities  for  7cm  chamber. 


Because  the  strength  of  the  magnetic  field  can  not  be  increased  as  chamber  size  is 
scaled  down,  the  Larmor  radius  increases  relative  to  the  size  of  the  chamber.  The 
primary  electrons  can  easily  cross  field  lines  reaching  all  regions  of  the  chamber.  For 
the  0.7cm  and  0.07cm  cases  then,  ionization  occurs  more  uniformly.  The  resulting 
plasma  density  distributions  are  shown  in  Figure  5-6  and  Figure  5-7.  A  lower  fraction 
of  ions  reach  the  beam  leading  to  decreased  performance  as  compared  to  the  7cm 
case. 
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axis 


Figure  5-6.  Plasma  densities  for  0.7cm  chamber. 


Figure  5-7.  Plasma  densities  for  0.07cm  chamber. 


For  the  cylindrical  geometry  analyzed  here  the  ratio  of  exit  area  to  total  chamber 
surface  area  is  0.1 1 .  It  is  clear  from  Table  5-1  that  as  chamber  size  is  scaled  down  the 
fraction  of  ions  reaching  the  beam  approaches  the  0.11  value.  This  indicates  that  the 
magnetic  field  does  not  significantly  confine  the  electrons  or  ions  at  the  smallest 
scales. 

5.3  SUMMARY  AND  CONCLUSIONS 

This  chapter  examined  possible  design  considerations  for  ion  microthrusters.  It  was 
shown  that,  although  desirable,  magnetic  fields  can  not  be  scaled  inversely 
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proportional  to  the  chamber  scale.  Even  the  best  practical  permanent  magnets 
provide  little  electron  or  ion  containment  below  chamber  sizes  of  ~1  mm.  At  those 
scales  the  fraction  of  ions  which  reach  the  beam  is  essentially  the  ratio  of  the  chamber 
exit  area  to  the  total  chamber  area.  The  next  chapter  applies  these  findings  to  the 
concept  of  a  linear  ion  microthruster. 
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Chapter  6 

Design  Considerations  for  the 
Linear  ion  Microthruster 


In  the  previous  chapter,  the  performance  of  cylindrical  ion  engines  as  they  are  scaled 
down  three  orders  of  magnitude  was  examined.  It  was  found  that  magnetic  field 
strengths  typical  of  permanent  magnets  had  no  effect  on  electron  or  ion  containment 
below  diameters  of  0.1mm.  This  resulted  in  the  fraction  of  ions  reaching  the  beam  and 
cathode  potential  surfaces  equaling  the  respective  surface  area  fractions.  This  result 
will  be  exploited  to  forego  consideration  of  magnetic  fields  in  analysis  of  the  linear  ion 
microthruster  concept. 


6.1  INTRODUCTION  TO  THE  LINEAR  ION  MICROTHRUSTER 

The  linear  ion  microthruster,  shown  in  Figure  6-1,  is  a  concept  for  a  micro-machined 
ion  propulsion  system  proposed  by  members  of  the  Jet  Propulsion  Laboratory  [1].  The 
ends  of  the  chamber  are  not  pictured  in  Figure  6-1  in  order  to  provide  an  unobstructed 
view  of  the  inside.  The  thruster  consists  of  several  linear  discharge  chambers  situated 
parallel  to  each  other.  Each  discharge  chamber  has  dimensions  of  approximately 
lOOjxm  x  300jim  x  10cm  although  these  sizes  are  only  conceptual.  The  walls 

separating  the  discharge  chambers  are  built  up  from  alternating  layers  of  conducting 
and  insulating  materials.  Contained  within  these  walls  is  the  ion  accelerator  system 
which  includes  the  ion  accelerator  grid  and  the  ion  decelerator  grid  separated  by 
layers  of  spray  coated  aluminum  oxide  of  sufficient  thickness  to  stand  off  the  high 
voltages  applied  between  the  layers.  Electrons  stripped  from  the  propellant  in  the 
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ionization  process  are  collected  by  the  anode  and  injected  into  the  ion  beam  by  an 
external  neutralizer.  The  linear  ion  microthruster  concept,  as  currently  conceived, 
does  not  contain  an  exit  grid.  Design  of  the  ion  acceleration  region  is  currently  on¬ 
going  at  the  Jet  Propulsion  Lab.  Consequently,  the  transparency  of  the  exit  region  to 
ions  and  neutrals  is  not  known  with  certainty.  In  this  chapter  then,  analysis  is  first 
performed  assuming  transparencies  typical  of  traditional  ion  engines  (i.e.  <J>,=0.86  and 

<j)o=0.26,  ref.  [2]).  Next,  the  analysis  is  repeated  with  ions  and  neutrals  treated  as 
though  they  can  pass  unhindered  through  the  exit  (i.e.  and  <1>0=1). 


ACCELERATOR  GRID 

SCREEN  GRID 
RE  “COIL" 

INSULATOR 

DISCHARGE  PLASMA  CHAMBER 


0.3  mm 


Figure  6-1.  Linear  Ion  Microthruster 

Brophy’s  model,  as  formulated  in  Chapter  4,  will  be  applied  in  the  analysis  of  this 
concept.  The  performance  of  an  ion  thruster  is  gauged  by  the  power  required  to 
produce,  but  not  accelerate,  an  ion  which  eventually  becomes  part  of  the  beam.  The 
average  energy  cost  to  produce  a  single  beam  ion,  %,  is  a  key  parameter  in  the 
measure  of  the  thruster’s  performance.  The  produce  a  more  efficient  thruster,  it  is 
necessary  to  reduce  the  value  of  eB  while  striving  to  maintain  a  high  propellant 
utilization. 
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6.2  PERFORMANCE  OF  LINEAR  GEOMETRY 

The  results  of  the  previous  chapter  demonstrated  that  as  size  is  scaled  down  several 
orders  of  magnitude  the  plasma  density  becomes  uniform.  Therefore,  the  problem  of 
finding  the  fraction  of  ions  reaching  the  beam  and  cathode  potential  surfaces,  fB  and  fc 
respectively,  is  a  simple  one  to  solve.  Their  values  correspond  to  the  surface  area 
ratios  of  the  beam  and  walls.  What  does  remain  to  be  found,  however,  is  the  primary 
electron  confinement  length,  le. 

6.2.1  Determination  of  Electron  Confinement  Length 

Arakawa’s  primary  electron  trajectory  simulator  was  modified  to  accommodate 
rectilinear  geometries.  Electron  trajectories  were  simulated  over  a  range  of  chamber 
height,  width,  and  lengths  to  determine  electron  confinement  lengths  for  each  case. 
Cathode  position  was  assumed  to  be  located  at  the  center  of  the  chamber  floor.  For  a 
given  chamber  width  to  height  ratio,  wyh,.,  confinement  length  was  found  to  be 
independent  of  chamber  length,  Lc.  Although  contrary  to  intuition,  the  reason  for  this 
can  be  seen  in  Figure  6-2  which  shows  an  electron  trajectory  for  a  linear  chamber.  If 
the  anode  is  located  in  the  sidewall  along  the  length  of  the  chamber,  then  the  path  of 
the  electron  (dotted  line)  before  it  reaches  the  anode  will  be  the  same  regardless  of 
chamber  length,  Lc. 
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The  anode  was  positioned  halfway  up  the  wall  with  a  constant  width  of  10pm.  All 
electrons  are  assumed  to  reflect  specularly  off  the  electrostatic  sheaths  surrounding 
the  chamber  surfaces  and  be  captured  upon  reaching  the  anode.  Arakawa’s  electron 
trajectory  simulation  showed  confinement  length  varying  with  chamber  width  to  height 
ratio  as  shown  in  Figure  6-3. 


Figure  6-3.  Variation  of  Electron  Confinement  Lengths  with  w^h,. 

From  Figure  6-3  it  can  be  seen  that  the  confinement  length  is  several  times  the 
chamber  dimensions  due  to  reflection  off  the  walls  prior  to  being  captured  by  the 
anode.  Further  the  confinement  length  decreases  as  chamber  width  to  height  ratio 
increases.  Because  the  anode  is  located  halfway  up  the  wall,  electrons  in  tall,  narrow 
(i.e.  low  Wj/hJ  chambers  tend  to  bounce  from  sidewall  to  sidewall  many  times  before 
reaching  the  anode.  The  optimal  confinement  length  of  2.1mm  appears  to  occur  at  a 
geometry  ratio  of  0.2.  This  corresponds  to  chamber  height  and  width  equal  to 
0.387mm  and  0.077mm,  respectively. 

6.2.2  Energy  Costs  per  Beam  Ion 

The  value  of  le,  found  in  the  previous  section  is  now  used  in  Brophy’s  model  to 
optimize  the  linear  microthruster  geometry  in  terms  of  height,  width,  length,  and 
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number  of  chambers.  Operating  parameters  were  assumed  to  be  consistent  with 
those  necessary  for  a  5  kg  microspacecraft  orbiting  as  part  of  a  distributed  satellite 
system.  In  Chapter  3,  it  was  found  that  the  maximum  tidal  accelerations  encountered 
by  such  a  “swarm”  of  spacecraft  would  be  on  the  order  of  0.001  m/s2.  Further  it  was 
shown  that  the  specific  impulses  required  are  typical  of  those  generated  by  ion 
engines. 


Figure  6-4  presents  the  methodology  used  to  analyze  the  performance  of  the  linear 
configuration.  Assuming  a  total  thrust,  Ftot,  of  5  mN  and  a  specific  impulse  of  3300s, 
the  total  mass  flow  rate  is  0.162  A-eq.  Knowing  the  chamber  width,  length  and  grid 
spacing,  the  allowable  mass  flow  rate  per  chamber  can  be  calculated  from  the  Child- 
Langmuir  1-D  space  charge  limit  (SCL) 


mion  = 


Ftot  _  ^ 

Ispg  9 


(6-1) 


where  Ac  and  Nc  are  the  chamber  exit  area  and  number  of  chambers,  respectively. 
Then  sgnd is  the  9aP  spacing  across  which  the  ion  beam  is  accelerated.  The  value  of 
Sgridforthe  linear  ion  microthruster  is  ambiguous  since  no  grid  physically  exists  at  the 
exit.  As  mentioned  earlier,  the  design  of  the  ion  acceleration  region  is  currently  under 
development.  Therefore,  for  this  analysis,  it  will  be  assumed  that  the  electric  potential 
which  accelerates  the  ions,  Vtot,  is  distributed  uniformly  between  the  side  walls  across 
a  gap  measuring  sgrid.  Allowable  mass  flow  per  chamber  can  be  divided  into  the  total 
mass  flow  for  the  thruster  yielding  the  required  number  of  chambers.  Values  of  Te  and 
£p  were  determined  through  an  iterative  process  based  on  the  energy  equilibrium  in 
the  discharge  chamber. 
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Figure  6-4.  Linear  analysis  methodology. 


Figure  6-5  shows  the  variation  of  beam  ion  production  cost  with  chamber  length  for  the 
case  where  transparencies  to  ions  and  neutrals  are  chosen  to  be  0.81  and  0.26, 
respectively.  The  ion  acceleration  region  was  limited  to  values  greater  than  10 
microns.  The  height  and  width  of  all  chambers  remained  constant  at  0.39  mm  and 
0.077  mm,  respectively.  At  a  chamber  length  of  0.70  mm,  15  chambers  were  required 
to  satisfy  space  charge  limitations.  As  length  increased,  fewer  chambers  were 
required  to  provide  the  necessary  mass  flow  rate.  Beam  ion  production  costs  tends  to 
decrease  as  the  number  chambers  decreases.  For  a  given  number  of  chambers, 
however,  beam  ion  production  cost  increases  as  length  increases.  This  is  due  to  an 
increase  in  exit  area  with  no  increase  in  mass  flow  rate  per  chamber.  The  trend  of  a 
slightly  decreasing  beam  ion  cost  reverses  dramatically  once  less  than  three 
chambers  are  needed.  For  that  situation,  beam  ion  cost  jumps  from  1000  to  5000  W/A 
as  the  length  of  the  chamber  increases.  For  relatively  short  chambers  (on  the  order  of 
chamber  height  and  width),  edge  effects  will  become  significant.  In  this  analysis, 
however,  no  attempt  was  made  to  model  the  effect. 


106 


The  beam  ion  production  costs  in  the  range  of  1100  to  5000  W/A  shown  in  Figure  6-5 
are  several  times  higher  than  typical  large  ion  engines.  This  is  due,  in  part,  to  a 
confinement  length  approximately  eight  times  chamber  dimensions.  This  value 
encompasses  the  ability  of  the  electron  to  reflect  off  walls  while  still  retaining  enough 
energy  to  ionize  neutrals.  However,  such  confinement  lengths  are  well  below  those  of 
large  ion  engines  (typically  20-40  times  chamber  dimensions) 
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Figure  6-5.  Relationship  between  beam  cost  and  number  of  chambers. 


Figure  6-6  presents  analysis  similar  to  that  just  performed  except  that  it  is  assumed  the 
chamber  exit  conditions  do  not  hinder  the  passage  of  ions  or  neutrals  (i.e.  <j>~1  and 

$0=1).  The  lower  neutral  density  in  the  chamber  results  in  a  jump  in  beam  ion 
production  cost.  Again  it  can  be  seen  that  beam  costs  decreases  slightly  as  the 
number  of  chambers  decreases.  Optimum  beam  ion  costs  occurs  at  2  chambers.  For 
a  given  number  of  chambers,  beam  ion  costs  increases  as  chamber  length  increases. 
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Figure  6-6.  Beam  ion  energy  cost  vs.  chamber  length 


6.2.3  Mapping  ion  Energy  Cost  to  Efficiency 

Before  the  performance  of  the  linear  microthruster  can  be  compared  to  other  systems 
in  terms  of  propulsion  system  power  and  mass,  it  is  necessary  to  put  the  cost  per 
beam  ion  value  in  terms  of  a  more  easily  used  parameter,  namely  thruster  efficiency,  tj. 

In  more  general  terms,  can  be  thought  of  as 

Total  Power -Useful  Power 

eB  = - — - 7 - = -  (6-2) 

JB 

where,  from  Chapter  4,  useful  power=JB(VD+Vnet)  and  total  power=Pin. 

Substituting  Pjn  from  equation  (6-2)  and  with  P^  =  JBVnet ,  the  efficiency,  r\,  of  the 
thruster  can  then  be  related  as: 

_  _  Pbeam 

P- 

rin 

TJ  =  - - - 

VD  +  Vnet  +eB 
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1 


(6-3) 


1  4  (VD  +gB)2e 
(goIsp)2mi 


Figure  6-7  shows  efficiencies  plotted  versus  specific  impulse  as  calculated  using 
equation  (6-3)  for  several  values  of  eb. 


Figure  6-7.  t|  vs.  Isp  at  Several  eB 


6.3  LINEAR  ION  MICROTHRUSTER  VS.  COLD  GAS  THRUSTER 

To  gauge  the  performance  of  the  linear  ion  microthruster,  a  comparison  is  made  with  a 
hypothetical  cold  gas  microthruster  for  a  range  of  thrust  times.  Using  equation  (6-3)  the 
lowest  Eg  in  Figure  6-6  corresponds  to  an  efficiency  q=16%.  The  cold  gas  thruster  is 

assumed  to  operate  at  a  Isp  =  75  s.  The  specific  power  of  the  ion  thruster  power  plant 
is  assumed  to  be  10  W/kg.  The  mass  of  both  the  cold  gas  and  linear  ion  microthrusters 
are  assumed  to  be  negligible  compared  with  the  mass  of  any  propellant  or  power  plant 
required. 

Figure  6-8  shows  how  propulsion  system  mass  varies  with  thrusting  time.  For  short 
thrusting  times  ,  the  mass  is  due  almost  exclusively  to  its  large  power  system.  The 
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efficiency  of  the  linear  ion  microthruster  is  so  poor  that  just  the  power  plant  makes  up 
approximately  55%  of  the  spacecraft  mass.  Obviously  this  is  not  acceptable.  For 
comparison,  a  more  efficient  engine  is  also  plotted.  Notice  that  at  5  years  thrust  time, 
the  more  efficient  ion  engine  reaches  the  levels  of  approximately  10%  propellant  mass 
and  10%  power  plant  mass.  This  is  consistent  with  the  findings  in  Chapter  3  that  the 
propulsion  system  of  spacecraft  in  clusters  should  have  lsp=3500  and  ri=65%.  The 

cold  gas  system,  whose  mass  is  made  up  entirely  of  propellant,  requires  less  mass 
than  the  linear  ion  engine  for  thrust  times  less  than  2  months.  As  thrust  times  increase, 
however,  the  lower  specific  impulse  of  the  cold  gas  system  forces  its  propellant  mass 
to  increase  quickly,  eventually  exceeding  the  mass  of  linear  ion  system.  A  linear  ion 
thruster  as  presently  conceived  is  not  appropriate  for  clusters  of  microsatellites. 


Figure  6-8.  Mission  Propulsion  mass  (excluding  thruster  itself) 

6.4  SUMMARY  AND  CONCLUSIONS 

The  performance  of  the  linear  ion  microthruster  was  examined  using  Brophy’s  model 
to  predict  energy  costs  per  beam  ion.  This  energy  cost  was  then  translated  into 
thruster  efficiencies.  For  linear  ion  microthrusters  on  the  scale  of  0. 1  mm  to  1  mm, 
magnetic  fields  can  not  be  used  to  increase  the  primary  electron  confinement  length. 
Path  lengths  are,  therefore,  limited  to  distances  on  the  order  of  chamber  dimensions. 


110 


Application  of  these  findings  to  linear  ion  microthruster  design  indicate  that  beam  ion 
production  costs  tends  to  decrease  as  the  number  of  chambers  decreases  for  a  given 
thrust  level.  Further,  for  a  given  number  of  chambers  beam  ion  cost  increases  as 
length  increases.  Shorter  chamber  lengths  are  more  efficient  due  to  a  lower  ratio  of 
wall  surface  area  to  exit  area.  This  results  in  a  lower  fraction  of  ions  neutralized  after 
contact  with  cathode  potential  surfaces.  For  short  thrusting  times  during  a 
microsatellite  cluster  mission,  a  cold  gas  system  is  feasible.  For  a  period  longer  than  a 
few  months,  the  mass  of  propellant  required  eliminates  a  cold  gas  system  from 
consideration.  Unfortunately,  the  linear  ion  microthruster,  as  currently  conceived,  is 
not  efficient  enough  to  perform  the  mission  for  any  length  of  time. 
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Chapter  7 

Summary  and  Conclusions 


The  previous  chapter  examined  design  considerations  relating  to  a  linear  ion 
microthruster  concept.  That  was  the  last  of  a  progression  of  topics  discussed  in 
this  thesis.  Those  topics  have  been  cost  modeling  of  distributed  satellite 
systems,  determination  of  propulsion  system  requirements  for  satellite  clusters 
and  swarms,  and  analysis  of  ion  micropropulsion  systems  for  use  in  swarms  of 
microsatellites.  This  chapter  summarizes  the  finding  of  this  thesis  and  reiterates 
the  conclusions  reached. 

In  Chapter  2,  the  total  costs  over  a  10  year  mission  life  were  calculated  for  three 
configurations  of  the  NPOESS  mission.  Distributing  the  primary  instruments 
among  three  smaller  satellites  significantly  increases  the  reliability  of  each 
individual  satellite.  This  increased  reliability  substantially  reduced  the  number 
of  ground  and  on-orbit  spares  required,  thus  reducing  mission  costs  compared 
to  a  single  large  satellite  configuration. 

From  this  analysis,  it  can  be  concluded  that  the  distribution  of  sensors  among 
several  satellites  is  appropriate  when  the  savings  from  increased  reliability 
outweigh  the  increased  cost  to  deploy  the  distributed  system.  Further, 
distribution  is  less  advantageous  for  systems  with  high  sensor  reliability. 

A  distributed  architecture  makes  sense  if  it  can  offer  reduced  cost  or  improved 
performance.  Because  the  performance  requirements,  and  the  associated 
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probability  of  satisfying  them,  are  embedded  in  its  calculation,  lifetime  cost  is  a 
useful  metric  for  architecture  analysis.  The  appropriateness  of  distribution  is 
highly  dependent  on  the  characteristics  of  a  specific  mission.  As  a  means  to 
reduce  mission  costs,  it  is  recommended  that  satellite  manufacturers  trade-off 
the  following  advantages  and  disadvantages  of  distribution  during  the  initial 
design  phase: 

•  Improved  resolution  corresponding  to  the  large  baselines  of  sparse 
apertures  formed  by  clusters  or  swarms  of  satellites. 

•  Lower  failure  compensation  costs  due  to  the  separation  of  important 
system  components  among  many  satellites. 

•  An  increase  in  system  complexity  leading  to  longer  development 
times. 

The  fundamental  conclusion  is  that  large  constellations  of  hundreds  or 
thousands  of  small  and  microsatellites  could  feasibly  perform  missions  currently 
being  carried  out  by  traditional  satellites  today. 

Chapter  3  examined  the  propulsive  requirements  necessary  to  maintain  the 
relative  positions  of  satellites  orbiting  in  a  local  cluster.  Formation  of  these  large 
baseline  arrays  could  allow  high  resolution  imaging  of  terrestrial  or 
astronomical  targets  using  techniques  similar  to  those  used  for  decades  in  radio 
interferometry.  A  key  factor  in  the  image  quality  is  the  relative  positions  of  the 
individual  apertures  in  the  sparse  array.  The  relative  positions  of  satellites  in  a 
cluster  are  altered  by  “tidal”  accelerations  which  are  a  function  of  the  cluster 
baseline  and  orbit  altitude. 

From  these  findings,  it  can  be  concluded  that  near  continuous  thrusting  by  a 
propulsive  system  is  necessary  to  maintain  the  relative  positions  of  the  satellites 
within  allowable  tolerances.  Satellite  mass,  volume,  and  power  constraints  limit 
reasonable  cluster  baselines  to  approximately  30  m,  300  m,  and  5000  m  at 
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1000  km,  10,000  km,  and  GEO  altitudes  respectively.  To  maintain  these  cluster 
baselines,  the  propulsive  system  must  operate  at  minimum  1^  and  efficiency 
equal  to  approximately  3000s  and  65%  respectively  with  a  thrust/spacecraft 
mass  ratio  of  approximately  1 5jiN/kg.  These  parameters  are  consistent  with 
those  of  ion  engines.  The  rest  of  this  thesis  examined  the  performance  of  micro 
ion  engines  with  the  purpose  if  using  them  to  maintain  the  relative  positions  of 
microsatellites  in  a  cluster  or  swarm. 

Chapter  4  formulates  an  analytical  model,  known  as  Brophy’s  Theory,  to  predict 
the  performance  of  cylindrical  ring-cusped  ion  engines.  The  model  is  used  to 
predict  the  performance  of  ring-cusped  ion  engines.  It  was  shown  that  energy 
cost  per  beam  ion,  eB,  can  be  formulated  as  a  simple  algebraic  equation.  Most 
of  the  parameters  in  Brophy’s  model  are  known  a  priori  given  the  choice  of 
propellant  gas,  chamber  geometry,  and  operating  conditions.  The  remaining 
parameters  can  be  determined  using  Arakawa’s  algorithm,  the  fundamentals  of 
which  were  also  described  in  this  chapter. 

Chapter  5  uses  Brophy’s  model  and  Arakawa’s  algorithm  to  examine  how 
performance  of  cylindrical  ion  engines  changes  as  their  scale  is  reduced 
several  orders  of  magnitude.  The  code  is  discussed  first  with  an  emphasis  on 
the  formulations  used  to  simulate  particle  behavior  within  the  chamber.  Next, 
an  argument  is  outlined  as  to  how  thruster  parameters  such  as  mass  flow  rate 
and  discharge  voltage  should  change  as  thruster  size  is  scaled.  The  effect  on 
performance  of  scaled  cylindrical  ion  thrusters  is  then  presented. 

Based  on  these  findings,  it  can  be  concluded  that  magnetic  fields  can  not  be 
scaled  inversely  proportional  to  the  chamber  scale.  At  those  scales  the  fraction 
of  ions  which  reach  the  beam  is  essentially  the  ratio  of  the  chamber  exit  area  to 
the  total  chamber  area.  Application  of  these  findings  to  linear  ion  microthruster 
design  indicate  that  beam  ion  production  costs  tends  to  decrease  as  the 
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number  chambers  decreases  for  a  given  thrust  level.  Further,  for  a  given 
number  of  chambers  beam  ion  cost  increases  as  length  increases. 

In  chapter  6,  the  performance  of  the  linear  ion  microthruster  was  examined 
using  Brophy’s  model  to  predict  energy  costs  per  beam  ion.  This  energy  cost 
was  then  translated  into  thruster  efficiencies.  For  linear  ion  microthrusters  on 
the  scale  of  0.1  mm  to  1  mm,  magnetic  fields  can  not  be  used  to  increase  the 
primary  electron  confinement  length.  Path  lengths  are,  therefore,  limited  to 
distances  on  the  order  of  chamber  dimensions. 

Application  of  these  findings  to  linear  ion  microthruster  design  indicate  that 
beam  ion  production  costs  tends  to  decrease  as  the  number  chambers 
decreases  for  a  given  thrust  level.  Further,  for  a  given  number  of  chambers 
beam  ion  cost  increases  as  length  increases.  Shorter  chamber  lengths  are 
more  efficient  due  to  a  lower  ratio  of  wall  surface  area  to  exit  area.  This  results 
in  a  lower  fraction  of  ions  neutralized  after  contact  with  cathode  potential 
surfaces.  Overall,  however  the  performance  of  the  linear  ion  microthruster  is 
very  poor.  For  short  thrusting  times  during  a  cluster  mission,  a  cold  gas  system 
is  feasible.  For  a  period  longer  than  a  few  months,  the  mass  of  propellant 
required  eliminates  a  cold  gas  system  from  consideration.  Unfortunately,  the 
linear  ion  microthruster,  as  currently  conceived,  is  not  efficient  enough  to 
perform  the  mission  for  any  length  of  time. 
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Abstract 

A  Particle-In-Cell  (PIC)  method  is  developed  and  applied  to  simulate  the  electron 
current  collection  by  a  positively  charged  tether  in  a  quiescent  unmagnetized  plasma 
under  the  Maxwellian  collisionless  condition.  We  compare  our  result  with  the  exact 
solution  to  validate  our  code.  This  simulation  is  performed  with  the  help  of  a  non- 
rectangular  grid  and  a  new  treatment  of  the  outside  boundary  condition.  The  error 
induced  by  a  non- rectangular  grid  is  calculated  and  the  effect  of  it  is  considered. 
The  outside  boundary  treatment  improves  the  accuracy  of  the  amount  of  current 
collected  by  the  tether.  A  very  small  ion  mass  is  used  and  it  is  verified  to  speed  up 
the  computation  considerably  without  loss  of  quality  in  the  result.  The  comparison 
with  the  exact  solution  shows  that  our  code  provides  good  qualitative  and  quantitative 
approximations. 
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Chapter  1 


Introduction 


In  tethered  satellite  technology,  it  is  important  to  estimate  how  many  electrons  a 
spacecraft  can  collect  from  its  surrounding  plasma  by  its  positively  charged  tether. 
The  analysis  is,  however,  very  difficult  because  of  the  small  but  significant  Geomag¬ 
netic  field  and  the  spacecraft’s  relative  motion  to  both  ions  and  electrons  [5].  One 
of  the  approaches  for  the  solution  to  this  problem  is  the  numerical  method.  In  the 
numerical  analysis  of  space  plasma,  one  of  the  most  reliable  methods  has  been  the 
Particle-In-Cell  (PIC)  method.  In  this  thesis,  we  develop  a  PIC  code  for  a  two  di¬ 
mensional  collisionless  plasma  without  magnetic  field. 

The  original  Particle-In-Cell  code  was  established  by  Birdsall  at  U.C.Berkeley.  [1] 
Using  a  rectangular  grid  and  finite-size  particles,  he  has  studied  the  effects  of  grid  size 
and  timestep  on  the  simulation  and  relevant  numerical  instabilities.  As  the  “bible” 
to  PIC  users,  his  publications  give  us  excellent  criteria  for  its  numerical  stability  and 
reliability.  The  concept  of  PIC  in  our  code  is  mostly  from  this  “bible”. 

The  usage  of  a  rectangular  grid  and  a  finite-size  particle,  however,  poses  some 
problems.  In  the  usual  numerical  application,  very  complicated  configurations  require 
non-rectangular  body  fitted  grids.  Being  defined  by  the  grid  size,  the  finite-size 
particle  can  not  maintain  its  constant  size  any  more.  In  this  thesis,  we  estimate  the 
error  induced  by  the  use  of  a  non-rectangular  grid  and  consider  the  effect  on  our 
simulation. 

In  order  to  validate  our  code,  we  apply  it  to  the  current  collection  by  a  cylindrical 
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tether  in  a  collisionless  unmagnetized  plasma,  near  the  boundary  of  Orbital  Motion 
Limited  (OML)  regime,  for  which  an  exact  solution  exists  [4].  The  OML  theory, 
which  will  be  further  discussed  in  the  next  chapter,  applies  in  the  limit  of  large  ratios 
of  Debye  length  to  radius. 

In  the  vicinity  of  the  tether,  there  is  a  region  called  “sheath”,  where  quasi¬ 
neutrality  does  not  apply.  That  is  to  say,  the  densities  of  ions  and  electrons  differ 
from  each  other  considerably.  In  order  to  reproduce  this  region,  we  follow  the  mo¬ 
tion  of  both  ion  and  electron  particles.  Since  the  mass  ratio  of  an  ion  to  an  electron 
is  very  large,  we  would  need  a  great  number  of  iterations  till  both  species  come  to 
have  converged  distributions,  without  having  electrons  travel  a  large  distance  in  one 
iteration,  which  would  induce  a  large  error  in  the  energy  conservation.  However,  the 
fact  that,  at  tether  potentials  much  greater  than  the  ion  temperature,  we  can  assume 
that  no  ion  is  absorbed  by  the  tether,  solves  this  problem  of  the  computational  cost. 
In  theory,  the  ion  density  at  an  arbitrary  point  does  not  depend  on  its  mass  when  its 
distribution  is  Maxwellian. 

Since  some  electrons  are  absorbed  on  the  surface  of  the  tether,  there  is  also  a 
region  called  “pre-sheath”  outside  the  sheath,  where  quasi-neutrality  prevails  but  the 
electric  potential  is  not  the  same  as  that  of  the  ambient  plasma  at  infinity. 

In  a  computation,  we  have  to  use  a  finite  region  to  calculate  field  quantities.  And, 
because  of  the  limited  memory  on  a  computer,  we  can  not  use  an  infinitely  large  grid 
to  include  the  pre-sheath  region.  Therefore  we  have  to  clip  a  computational  region 
out  of  this  infinitely  large  space,  and  determine  the  outer  boundary  conditions  by 
considering  only  the  quasi-neutrality  and  the  collisionless  nature  of  the  plasma.  The 
limitation  of  the  number  of  particles  available  in  a  computer  gives  fluctuating  local 
boundary  conditions.  To  avoid  this,  we  use  a  spatially  averaged  boundary  condition. 

The  goal  of  this  thesis  is  to  establish  a  code  to  simulate  a  collisionless  unmagne¬ 
tized  plasma  in  and  near  the  sheath  region.  Based  on  this  code,  we  plan  to  include 
later  all  the  phenomena  encountered  by  a  tethered  satellite,  such  as  Geomagnetic 
field  and  plasma  cross-flow. 
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Chapter  2 


Orbital  Motion  Limit  (OML) 
Current 


Current  collection  by  spherical  and  cylindrical  probes  (tethers)  was  first  analyzed 
by  Lagmuir  and  Mott-Smith  [6]  ,  who  named  the  thin  cylinder  limit,  ’Orbital  Mo¬ 
tion  Limit  (OML)’.  When  OML  theory  applies,  namely,  when  the  ratio  of  the  probe 
(tether)  radius  to  the  Debye  length  of  the  plasma  is  so  small  that  the  shielding  be¬ 
comes  unimportant,  the  number  of  electrons  absorbed  by  the  probe  is  determined 
from  energy  and  angular  momentum  considerations  alone. 

The  OML  limit  can  be  described  in  terms  of  the  effective  potential.  [3]  Let  J  and 
E  be  the  angular  momentum  and  the  energy  of  a  particle,  respectively.  From  the 
energy  conservation  and  the  angular  momentum  conservation  of  an  electron  in  two 
dimensions,  although  the  velocity  uj|  along  the  cylinder  axis  can  be  nonzero,  we  have 

E  =  ^me  ( v l  +  vf)  +  q<f>  (2.1) 

J  —  mervg  (2.2) 

where  r  is  the  distance  from  the  probe  center,  me  the  electron  mass,  q  the  electron 
charge  of  an  electron,  (f>  the  local  potential,  vr  the  radial  velocity  component,  and  v$ 
the  azimuthal  velocity  component.  Substituting  equation  (2.2)  into  equation  (2.1), 
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we  have 


v2r  =  —  [e  -  q<t>  - 

me  \  2 mer2) 


(2.3) 


In  order  for  a  particle  to  reach  the  surface  of  the  probe,  the  right-hand  side  of  equation 
(2.3)  must  be  positive  not  only  at  the  surface  of  the  probe,  but  also  all  along  the  path 
from  infinity  to  the  surface.  To  consider  the  particle  motion  from  the  one  dimensional 
viewpoint,  the  effective  potential  defined  by 


U  =  q<j>  + 


J2 

2  mer2 


(2.4) 


should  be  considered.  Substituting  the  effective  potential  (2.4)  into  (2.3),  we  have 


v2  =  ^-(E-U). 

7Tte 


(2.5) 


By  taking  the  effective  potential  as  a  normal  potential,  we  can  treat  the  2-dimensional 
particle  motion  as  the  1-dimensional  case.  Fig(2-1)  illustrates  two  limits  regarding  the 
effective  potential.  Assume  that  the  probe  is  on  the  left  of  the  figure.  When  the  sheath 
is  thin  (Langmuir  Limit),  the  second  term  of  equation  (2.4)  becomes  dominant  near 
the  probe  and  v(r)  has  an  intermediate  minimum  value.  For  some  attracted  particles, 
this  bump  in  the  effective  potential  prevents  them  from  reaching  the  surface  even  if 
they  have  enough  energy.  When  the  sheath  is  thick  (OML  limit),  the  first  term  in 
equation  (2.4)  becomes  dominant  throughout  the  region,  and  the  electric  potential  is 
large  enough  to  overwhelm  the  bump  in  the  effective  potential.  Therefore  the  effective 
potential  becomes  monotonous,  and  the  only  requirement  for  a  particle  to  reach  the 
surface  is  to  have  a  positive  value  of  the  right-hand  side  of  equation  (2.3)  at  the  probe 
surface. 

Electrons  absorbed  by  the  probe  should  be  accelerated  by  the  field  force  up  to  a 
certain  total  velocity  toward  the  probe.  Therefore,  in  terms  of  energy,  it  is  equivalent 
to  say, 

^me(v2  +  vg)  +  #p  >  0  (2.6) 

where  <f)p  is  the  probe  potential. 
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Figure  2-1:  The  effective  potential  for  the  thick  and  thin  sheath 


In  the  absence  of  collisions,  the  solution  to  Vlasov’s  equation  must  have  the 
Maxwell-Boltzman  form  for  any  velocity  that  does  occur: 


where  n <*,  is  the  density  at  infinity,  k  the  Boltzmann  constant,  the  temperature 
at  infinity  and  i/||  a  velocity  component  parallel  to  the  cylindrical  probe. 

At  the  surface  of  the  probe,  only  electrons  which  satisfy  the  equation  (2.6)  can 
exist  and  be  counted  for  the  current  collection.  The  current  density  into  the  probe  is 
given  as 


3  =  Q 


=  q 


f!L 
L 


dv\\[ 

_  r 

qriooCoc 


fedv 


rw/2 

/  v±fevj_cos9d9dv± 
J—tt/2 


2y/n 


y^l^pl/me  J-1 r/2 

m+>/ze^erfc(ML 

MU  2  MV  kTra 


(2.8) 

(2.9) 

(2.10) 


11 


g^ooCpo  l\q<t>p\ 

20F  V  kT^ 

since,  in  the  limiting  form,  x  -¥  oo,  we  have 


(2.11) 


erfc(. 


x)  =  1  ~  erf(x)  =  ~  [  e~t2dt 

V7T 


2  e“x2 


\/^  x 

where  C'oc  is  the  random  thermal  velocity  given  as 


8  kT« 


nmP 


Therefore,  when  »  1,  the  current  density  (2.13)  becomes 


(2.12) 

(2.13) 


(2.14) 


3  = 


qnco 

7T 


(2.15) 


which  is  independent  of  electron  temperature,  Te.  Note  that  equation  (2.10),  and 
hence  equation  (2.15)  are  independent  of  the  shape  of  the  cylinder’s  cross  section  (as 
long  as  OML  conditions  prevails). 

This  limiting  value  of  current  density  is  used  in  this  thesis  as  one  of  the  criteria  in 
the  validation  of  our  code.  To  see  whether  our  code  simulates  plasma  behavior  well 
in  the  vicinity  of  the  OML  regime,  we  examine  several  Debye  ratios,  rp/dDebye,  for 
the  case  that  the  ratio  of  the  probe  potential  to  the  electron  temperature,  yp  = 
is  25,  which  case  Laframboise  has  exactly  computed.  The  solutions  are  available  as 
analytical  fits,  which  we  quote  in  this  thesis.  [4,  2] 
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Chapter  3 


Numerical  Method 
(Particle-In-Cell) 

In  this  chapter,  the  Particle-In-Cell  (PIC)  method  used  here  is  explained  in  detail. 
First,  we  introduce  the  structure  and  mechanism  of  the  PIC  method.  Secondly,  we 
describe  the  model  for  our  simulation  and  the  governing  equations  in  a  nondimen- 
sionalized  form.  Thirdly,  we  discuss  the  problems  with  the  PIC  method,  which  have 
occured  in  our  application.  Finally,  we  show  the  results  from  this  simulation  and 
compare  them  with  the  exact  solutions. 

3.1  Particle-In-Cell  (PIC) 

The  Particle-In-Cell  (PIC)  method  has  been  very  successful  in  the  simulation  of  col¬ 
lisionless  plasmas.  In  PIC,  many  particles  are  distributed  in  phase  space.  That  is,  a 
particle’s  motion  is  described  by  its  position  and  velocity.  In  kinetic  theory,  this  par¬ 
ticle  distribution  is  defined  as  a  distribution  function  and  governed  by  the  Boltzmann 
equation.  The  Boltzmann  equation  with  no  collisional  term  on  its  right-hand  side  is 
given  as  follows.  (Vlasov’s  equation) 


In  an  actual  computation,  the  number  of  particles  available  is  much  less  than 
that  in  reality.  This  fact  requires  us  to  introduce  the  concept  of  a  “superparticle”, 
corresponding  to  a  group  of  real  particles.  One  superparticle  contains  many  real 
particles,  and  as  many  particles  as  another. 

To  describe  the  motion  of  the  superparticle,  we  need  to  know  velocity  and  the 
force  acting  on  it.  The  force  acting  on  a  superparticle  can  be  calculated  by  consid¬ 
ering  all  Lorentz  forces  caused  by  the  other  superparticles.  However,  this  calculation 
is  computationally  too  expensive.  Instead  of  doing  so,  PIC  uses  a  grid  on  which 
Maxwell’s  equations  are  solved  to  give  the  electric  field,  which  is  then  interpolated 
to  the  position  of  each  superparticle.  As  the  name  “Particle-In-Cell”  implies,  in  a 
computational  domain,  a  superparticle  moves  through  a  grid  or  a  cell,  regardless  the 
position  of  grid  nodes.  A  PIC  code  method  consists  of  four  processes  as  described 
below. 

At  each  time  step,  the  electric  charge  density  on  each  node  is  estimated  from 
the  positions  of  all  superparticles.  This  first  process  is  called  “charge  assignment”. 
Then  on  a  grid,  the  electric  potential  and  electric  field  are  computed.  We  use  a 
finite  difference  method  in  this  second  process;  especially  to  solve  Poisson’s  equation, 
we  use  Successive  Line  Over  Relaxation  (SLOR).  Poisson’s  equation  to  relate  electric 
potential  to  charge  density  is 

V20  =  (3.2) 

where  p  =  |g|(nj  —  ne)  is  the  electric  charge  density  ,  and  the  electric  field  is 

E  =  -V<£.  (3.3) 


if  ®  can  be  neglected. 

After  computing  E  on  a  grid,  the  electric  field  is  interpolated  onto  each  super¬ 
particle’s  position,  and  the  corresponding  force  and  acceleration  of  a  superparticle 
is  calculated.  This  third  process  is  called  “interpolation”.  The  first  and  third  pro¬ 
cess  involve  the  same  weighting  method  to  avoid  the  so-called  “self  force”.  Once  we 
know  the  acceleration,  a  leapfrog  method,  the  final  process,  updates  the  velocity  and 
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position  of  each  superparticle  as  follows; 


» 


x(«+l/2) 

new 


-  v(n 

—  Wold 


-i)  +  gE<"-V2) 


At 


m 


=  x(?71/2^  +  v(n)  At 

^old  '  v*j*>ifi*-*t'* 


(3.4) 

(3.5) 


This  completes  one  iteration  in  a  PIC  calculation.  One  cycle  of  a  PIC  is  shown 
skematically  in  fig  3-1. 


3.2  Simulation  Model 

This  section  explains  the  simulation  model.  We  wish  to  calculate  the  electron  cur¬ 
rent  collection  by  a  positively  charged  cylindrical  tether  in  a  quiescent  unmagnetized 
plasma  in  a  Maxwellian  collisionless  condition.  For  simplicity,  we  first  nondimension- 
alize  the  governing  equations. 

Non-dimensionalization 

Before  we  consider  the  non-dimensionalization,  we  should  know  what  equations  are 
involved  in  our  simulation.  As  we  showed  before,  the  leapfrog  method  uses  equations 
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(3.4)  and  (3.5).  To  solve  for  electric  field  from  electric  charge  density,  we  use  Poisson’s 
equation  (3.2)  and  Maxwell’s  equation  (3.3).  Essentially,  these  four  equations  are  the 
governing  equations.  Since  we  are  not  directly  solving  Boltzmann’s  equation  (3.1),  it 
is  not  considered  as  a  governing  equation.  The  distribution  function,  however,  should 
be  considered.  A  Maxwellian  distribution  function  is  used  to  calculate  the  number  of 
particles  replenished  into  the  computational  domain  at  each  timestep  and  the  density 
at  the  outside  boundary. 

In  non-dimensionalizing  the  governing  equations,  we  use  reference  values  as  fol¬ 
lows; 


Length 


.  /  _  j  _  UKT, 

•  bef  ~  “ Debye  ~  \f 


Time  :  tref  —  10/cup  (up  —  y  ) 

Potential  :  <j)ref  = 

Density  :  nref  =  71^/100 

Velocity  :  vref  =  lref/tref  =  vT/10  =  10 

Distribution  function  :  fref  =  nref/vfej. 

As  is  discussed  later,  somehow  strange  reference  values  seen  here  is  totally  due  to  the 


computational  limitations. 


Substituting  these  reference  values,  we  nondimensionalize  equations  (3.2),  (3.3), (3.4), (3.5) 
and  (2.7). 


V2«£  =  -(n i  -  ne)/ 100  (3.6) 

E  =  -V4>  (3.7) 

^ new  =  V0,d  ±  100(— )Edi  (3.8) 

m 

X new  X old.  ^new^t  (3.9) 

f  =  10(2 hffi exp  -  *)  <31°) 


where  a  hat  Q  indicates  a  nondimensional  quantity.  In  equation  (3.8),  the  minus 
sign  is  taken  when  the  particle  is  an  electron  (m  =  me),  and  the  plus  when  an  ion 

(m  =  nii). 

In  this  simulation,  we  use  the  parameters  in  Table  3.1.  In  deciding  these  values, 
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Electron  Temperature 

EHi 

Ion  Temperature 

fcfeMSl 

Electron  density 

he  =  100 

Ion  density 

hi  =  100 

Table  3.1:  Parameters 


we  first  consider  the  number  of  particles  available  in  a  computer  memory.  In  order 
to  include  the  sheath  region  completely,  we  have  the  radius  of  the  computational 
domain  as  1  bd^ebye-  To  be  consistent  with  the  governing  equations  and  these  non¬ 
physical  values,  we  calculate  other  nondimensionalized  variables,  starting  with  the 
Debye  length  equal  to  unity, 


d 


Debye 


(3.11) 


Since  the  area  of  the  computational  domain  is  ( 15d£>ebye)2n ,  we  need  approximately 
heti(15dDebye)2Tr  particles  for  each  species.  This  number  is  limited  by  the  computer. 
In  our  case,  the  number  of  particles  available  is  about  200, 000.  In  order  to  run  the 
simulation  with  as  many  particles  as  possible  but  less  than  this,  we  set  ne,i  =  100. 
From  equation  (3.11),  we  also  have  the  temperature,  T^e  =  100.  Consequently,  we 
have  the  thermal  velocity, 

V^00  (3.12) 


and  the  plasma  frequency, 


-  /<Z2n 

“P  =  X—tref 

V  em  J 

=  vToo. 


(3.13) 
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As  we  can  see  now,  the  factors  seen  in  the  reference  values  are  determined  by  the 
computational  limitations. 

We  should  note  that  we  use  hypothetically  light  ion’s  mass,  which  is  the  same  as 
an  electron’s  mass.  We  discuss  the  validity  of  this  hypothesis  and  its  effect  on  the 
computational  cost  later.  Physically,  since  almost  no  ions  are  lost  to  the  probe,  all 
velocities  are  possible  everywhere,  and  their  distribution  is,  in  fact,  a  full  Maxwell- 
Boltzmann  distribution,  so  that  the  ion  density  is  simply  n*  =  exp(-j^),  which 

does  not  depend  on  at  all. 

Mesh 

In  our  code,  the  Debye  length  is  first  determined.  Therefore,  for  various  cases  of 
Debye  ratio,  i.e.  the  ratio  of  the  tether  radius  to  the  Debye  length,  different  meshes 
are  used.  The  mesh  used  for  the  case  of  the  Debye  ratio  equal  to  1,  =  1,  is  shown 

in  Figure  3-2.  Each  mesh  size  in  the  radial  direction  is  kept  to  be  a  half  of  the  Debye 
length  and  in  the  azimuthal  direction  mesh  size  is  kept  less  than  the  Debye  length, 
which  should  avoid  numerical  instabilities. 

The  main  purpose  of  this  simulation  is  to  calculate  how  many  particles,  mainly 
electrons,  are  collected  by  a  cylindrical  tether  in  a  Maxwellian  collisionless  condition. 
Particles  are  counted  and  absorbed  when  they  reach  the  surface  of  the  tether.  From 
outside  of  the  computational  domain,  where  plasma  is  assumed  to  be  Maxwellian, 
electrons  and  ions  are  replenished  into  the  domain  with  velocity  and  position  calcu¬ 
lated  from  the  Maxwellian  distribution  function  using  a  random  number  generator. 


3.3  Problems  in  PIC 

As  we  have  applied  the  PIC  method,  we  have  been  confronted  with  some  problems. 
In  this  section,  we  discuss  those  problems  and  some  solutions  to  them. 
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3.3.1  Deformed  Grid 


As  we  consider  some  practical  configurations  in  engineering  applications,  it  is  almost 
impossible  to  do  it  with  only  a  rectangular  grid.  The  more  complicated  the  configura¬ 
tion  becomes,  the  more  necessary  it  is  to  use  a  deformed  body-fitted  grid.  Concerning 
this  problem,  we  can  think  of  two  factors  of  a  deformed  grid,  which  are  likely  to  cause 
problems,  namely  the  area  of  the  cell  and  the  degree  of  deformity  from  a  rectangle. 

First,  we  describe  the  mechanism  of  the  “assignment”  and  “interpolation”  pro¬ 
cesses.  In  the  original  PIC  method  based  on  a  rectangular  grid,  electric  charge  density 
is  first  calculated  from  the  particle  position  and  assigned  to  four  nodes  according  to 
the  area-weighting  function.  The  area-weighting  function  is  defined  as  follows. 


Uniform  rectangular  grid 

A  particle  with  electric  charge  q  is  located  at  (x,  y)  in  a  cell  of  area  A,  which  is  defined 
by  four  nodes  fa,  yt),  (re*  +  Ax ,  yi),  (xj,  yi  +  Ay)  and  ( Xi  +  Ax ,  yi  +  Ay).  That  is,  the 
area  A  is  given  as 

A  =  Ax  Ay.  (3.14) 

First,  we  calculate  the  electric  charge  density,  by  dividing  q  by  the  area,  Ax  Ay. 


P  = 


Q 

AxAy 


(3.15) 


Next,  we  split  this  charge  density  into  four  segments,  which  are  proportional  to 
the  area  demarcated  by  lines  parallel  to  grid  edges.  Thus,  to  point  A(xi,yi),  charge 
density 


=  VPFCG  (si+1  -  x)(yi+ 1  -  y) 

PA  PQABCD  9  Ax2  Ay2 

is  assigned.  Likewise,  to  the  other  points,  we  assign  charge  density  as 


(3.16) 


UPGDH_  _  (x  -  Xi)(yi+i  -  y) 
P  D  ABC  D  9  Ax2  Ay2 

OPHAE  _  (x  -  Xj)(y  -  yi ) 
p  □. ABC D  ~  9  Ax2 Ay2 


(3.17) 

(3.18) 
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Figure  3-3:  A  uniform  rectangular  grid 


UPEBF  _  (xi+i  -  x)(y  -  yi ) 
PUABCD  ~  q  Ax2  Ay2 


(3.19) 


This  assignment  process  is  applied  to  all  particles  in  the  cell  ABCD. 

We  next  estimate  the  order  of  accuracy  of  this  assignment  method.  To  exclude 
the  effects  of  an  insufficient  number  of  particles,  we  assume  that  there  is  a  sufficient 
number  of  particles  in  a  cell.  This  condition  is  expressed  by  requiring  the  density  to 
be  a  continuous  function.  When  we  consider  a  small  element  in  the  cell,  the  condition 
requires  that  there  be  still  many  particles  in  it. 

Let  f(x,y)  be  the  superparticle  density  function  at  (x,y),  and  dxdy  be  the  area 
of  the  small  element.  The  condition  to  have  a  sufficient  number  of  particles  requires 


f(x,  y)dxdy  >  1 


(3.20) 


Assuming  a  large  enough  number  of  particles  in  a  cell  and  using  the  Taylor  series 
expansion,  we  estimate  the  order  of  accuracy  at  the  node  A.  One  assignment  process 
is  applied  at  once  to  all  particles  in  the  small  element.  Since  the  element  is  taken 
to  be  very  small,  we  can  consider  that  there  are  f(x,y)dxdy  particles  at  (x,  y).  The 
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D  C 


A(Xi,yi)  X  B 


Figure  3-4:  A  small  element  in  a  cell 


assigned  charge  density  at  A  from  those  particles’  point  ( x ,  y)  is  given  by 


dPA(x,y)  =  -f(x,y)dxdy. 


(3.21) 


Assigning  all  particles  in  a  whole  cell  to  the  node  A,  we  have, 


rxi+i  rm+i 

Pa  —  /  dpA{x,y) 

Jxi  Jyi 

[Xi+1  fVi+i  gfoi+i  -s)(&+i  -y) 

4  4  Ax2Ay2  H  ’* 

=  f1  [lq(l-X)(l-Y)g(X,Y)dXdY 
Jo  Jo 


f(x,  y)dxdy 


(3.22) 

(3.23) 

(3.24) 


where  X  =  Y  =  ^  and  g(X,Y )  =  f(xt  +  AxX,yi  +  AyY).  Expanding 
g(X ,  Y)  around  the  point  (0, 0), 


g(X,Y)  =  </(0,0)  +  ||  X+^  K  +  - 

(0,0)  (0,0) 

=  f{xi,Vi)+y~  AxX+^~  AyY 
dx  (*,.»)  ^  (««*) 


(3.25) 


(3.26) 
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and  substituting  this  into  equation  (3.24),  we  have 


PA  =  q  fl  hnx^Xl-XXl-Y) 
J  0  J  0 


+ 


dg_ 

dX 


X(l-X)(l-F)  + 


(0,0) 


dg 

dY 


Y(  1  -  X)(l  -  Y)  +  •  •  ytXdY{ 3.27) 

(0,0) 


-  lftx.v.)+±M. 

AJ{xuyt)+  l2dx 


A  ,  Q  df 

AX+nTy 


A y  + 

(Xi,yt) 


(3.28) 


Performing  the  same  calculation  for  all  particles  in  the  other  cells  which  surround  the 
node  A,  we  can  cancel  the  2nd  and  3rd  terms  and  get 


Pa  =  qf(xi ,  Vi)  +  0( Ax2,  Ay2). 


(3.29) 


This  shows  that,  as  long  as  there  are  a  sufficient  number  of  particles,  this  area¬ 
weighting  assignment  method  provides  2nd  order  of  accuracy. 

Rectangular  grid  with  different  cell  sizes 

Next,  following  the  same  procedure,  we  consider  the  effect  of  rectangular  cells  of 
different  sizes  on  the  assignment  method  (fig  3-5).  When  the  cells  are  of  different 
sizes  but  still  rectangular,  the  2nd  and  3rd  terms  in  equation  (3.28)  do  not  cancel  out 
after  the  summation  of  corresponding  terms  from  other  cells.  Instead  of  the  second 
order  of  accuracy,  we  get 


Pa  =  gfixuVi) 


Q  df_ 

6  dx 


( Aii 

+0( Ax2,  Ay2), 


A*2)  +  f  f 

6  dy 


(Ayi  -  A y2) 


(3.30) 


which  provides  1st  order  of  accuracy  only. 

As  we  can  see  from  the  2nd  and  3rd  terms  in  equation  (3.30),  using  almost  the 
same  shape  and  size  as  those  of  adjacent  cells,  we  can  make  these  terms  quite  small, 
and  make  this  method  closer  to  second  order  of  accuracy. 
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Deformed  grid 

Next  we  consider  the  deformity  of  a  grid.  Since  we  use  a  linear  interpolation  [8,  7] 
from  a  deformed  grid  to  a  square  grid  on  which  the  area- weighting  is  performed, 
we  can  not  avoid  a  problem  due  to  this  deforming.  In  our  code,  we  use  the  typical 
tri-linear  interpolation  from  the  physical  domain  ( x ,  y)  to  the  computational  domain 
(r,  s)  given  as 


(3.31) 

(3.32) 

where 

Xq  Xj  j 

Vo  =  Vi,j 

-^1  *t+lj  Xij 


x  =  x0  +  x\r  +  x2s  +  x$rs 
V  =  yo  +  yir  +  y2s  +  y3rs 
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Figure  3-6:  Transformed  square  grid  (left)  and  the  original  deformed  grid  (right) 


2/i  —  Vi+l,j  Ui,j 

X2  =  Xij+i  Xij 

2/2  =  UiJ+1  ~~  Vi,j 

X3  =  ^i+l  j'+l  +  Xij  Xi+l,j  -^jll 

2/3  =  Vi+hj+l  +  Vi,j  -  Vi+l,j  -  Vij+l 

The  effect  of  the  linear  interpolation  is  understood  when  we  look  at  Figure  3-6. 

In  a  deformed  cell,  the  grid  density  is  sparser  at  the  top  and  denser  at  the  bottom. 
But  its  corresponding  square  cell  has  a  completely  uniform  density  distribution.  This 
means  that  if  we  start  for  example  with  a  particle  distribution  which  is  uniform  in  the 
grid  on  the  left,  it  will  not  be  uniform  in  the  one  on  the  right,  and  vice  versa.  Although 
we  have  no  particular  treatment  done  to  solve  this  problem  in  our  simulation,  as  seen 
in  the  grid  we  use  (fig.  3-2)  ,  each  cell  is  almost  rectangular.  Therefore,  no  serious 
error  can  be  considered  to  be  incurred  due  to  the  deformity  of  the  cell.  Moreover,  in 
the  grid  no  pair  of  adjacent  cells  has  a  significant  difference  in  size  and  shape.  This 
makes  the  assignment  method  to  be  close  to  2nd  order  accurate. 
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3.3.2  Boundary  Condition 

In  this  section,  we  consider  the  outside  boundary  condition  in  our  simulation.  The 
importance  of  this  boundary  condition  stems  from  the  fact  that  some  electrons  are 
absorbed  by  a  tether.  We  assume  here  that  no  ion  is  absorbed  by  the  tether,  because  of 
its  very  high  positive  potential.  Figure  3-7  illustrates  the  overall  flow  of  electrons  and 
ions.  Due  to  this  partial  absorption  of  electrons  by  the  tether,  the  electric  potential 
at  the  computational  outside  boundary  can  not  be  zero  with  respect  to  the  ambient 
plasma.  If  it  were  zero,  the  electron  density  would  be  less  than  the  ion  density,  and 
thus  it  would  violate  the  quasi-neutrality  outside  the  sheath.  To  maintain  the  quasi¬ 
neutrality,  the  electron  potential  at  the  outside  boundary  should  be  more  than  zero. 
Positive  potential  attracts  more  electrons  and  fewer  ions.  This  boundary  condition  is 
formulated  as  follows. 

Let  <f)  be  the  electric  potential  at  an  arbitrary  point  including  the  outside  boundary. 
Assuming  that  ions  are  singly  charged,  we  have  the  quasi-neutrality  equation  as 

| ne  —  rii\  ne.  (3.33) 


In  the  computation,  we  use  this  condition  in  the  form  of 

ne  =  j'ij.  (3.34) 

However  this  does  not  allow  us  to  transform  Poisson’s  equation  to  Laplace’s  equation 
by  equating  the  source  term  to  zero,  because  the  small  difference  \q\(ni—ne)  is  divided 
by  the  small  quantity  e0,  leaving  V2(f>  indeterminate.  As  the  plasma  approximation 
claims,  plasma  tends  to  neutralize  itself  by  imposing  ne  =  n*.  Therefore  we  impose 
the  condition  (3.34)  on  the  outside  boundary,  and  solve  Poisson’s  equation  inside  that 
boundary  only  with  non-zero  source  term  on  its  right-hand  side. 

To  determine  the  boundary  condition,  we  need  the  potential,  which  is  calculated 
as  follows.  The  assumption  that  no  ions  are  absorbed  by  the  tether  enables  us  to  cal¬ 
culate  the  ion  density  at  any  point.  Given  the  ion  temperature,  7*  and  the  potential, 
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4>,  and  integrating  the  Maxwell-Boltzmann  distribution  function  in  velocity  space,  we 
have  the  ion  density  as 

rii  =  rioc  exp(-® )  (3.35) 

On  the  other  hand,  we  can  calculate  only  the  density  of  inbound  electrons1.  Since 
some  electrons  are  absorbed,  we  do  not  know  the  limits  of  outbound  electron  distri¬ 
bution  function  in  phase  space.  We  denote  the  density  of  those  outbound  electrons 
as  nf*.  As  for  the  inbound  electron  density,  we  can  calculate  it  by  integrating  the 
Maxwellian  distribution,  since  all  those  electrons  can  be  tracked  back  to  infinity, 
where  the  Maxwellian  distribution  prevails.  The  inbound  electron  density  is  given  at 
any  radius  by 


n: 


n7r/2  roc 

72  J\ 


Tt/2J  \/\q\<l)lm{ 


+  tj)  +  q£ 
kT 


^oo 

~2~' 


v cos  OdvdOdvz 
(3.36) 


Substituting  (3.35)  and  (3.36)  into  (3.34),  Then,  we  have 

„r«  +  !|2  =  „0oeXp(-i|M)  (3.37) 

Given  the  outbound  electron  density  at  the  outside  boundary,  equation  (3.37)  allows 
us  to  calculate  the  potential  there. 

The  outbound  electron  density,  nf1*  is  calculated  computationally  by  considering 
the  outbound  flux,  T.  The  flux  through  the  boundary  is  given  by 


r  =  (3.38) 


where  is  the  flow  velocity  due  to  the  outbound  electrons,  that  is,  the  average 


1In  this  thesis,  we  call  a  particle  which  is  coming  into  the  computational  domain,  an  “inbound” 
particle,  and  one  going  out  of  the  domain,  an  “outbound”  particle. 
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Ion  Flow 


Figure  3-7:  Electron  and  Ion  Flow 


velocity  normal  to  the  boundary,  given  as 


v, 


out 


Eil  W ei  •  ru 
k 


(3.39) 


where  k  is  the  number  of  electrons  counted  as  they  cross  the  boundary,  n,  the  normal 
vector  to  the  boundary  and  wei  the  particle  velocity.  The  number  of  electrons  which 
go  out  of  the  domain,  k,  is 

k  =  TdtS  (3.40) 


where  dt  is  the  timestep  and  S  is  the  area  of  the  outside  boundary.  Substituting 
(3.40)  and  (3.39)  into  (3.38),  we  have  the  outbound  electron  density  as 


n°f  = 


k2 


Ei=i 


(3.41) 


Now  we  are  ready  to  calculate  how  many  electrons  and  ions  are  to  be  replenished 
at  each  timestep.  The  number  of  those  particles  is  calculated  by  multiplying  the  flux 
by  the  timestep  and  the  area  of  the  outside  boundary.  The  number  of  electrons  to  be 
replenished,  ke,  is  given 


rco  riT/2roo  /  m  \3/2  /  ime(v 2  +  V2)  +  q(f)\  2  ... 

ke  =  I  /  /  , _ rioo  (  — —  )  exp  (  — - - — -  v2  cos  9dvd0dvzSdt 

J—  <xJir/2J  y/\q\ij>/m^  X’l'ItkT )  \  kT  J 


TIqqCqo  \  2  j 

4  1  \/tt  V  kT 


+  exp 


(*)"#} 


Sdt 


(3.42) 


where  c0 Q  is  the  random  thermal  velocity  given  in  equation  (2.14).  It  should  be  noted 
that  this  is  the  same  as  equation  (2.10),  except  for  <f>  instead  of  <f>p.  And  the  number 
of  ions  to  be  replenished,  ki,  is  given 


ki 


v2  cos  9dvd9dvzSdt 


(3.43) 


In  the  simulation,  a  random  number  generator  is  used  to  locate  the  place  to  be 
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replenished  and  provide  the  particle  with  a  velocity  estimated  from  the  Maxwellian 
distribution  function. 


3.3.3  Computational  Cost 

In  order  to  resolve  the  non-neutral  domain,  namely  a  sheath  domain,  both  electrons 
and  ions  should  be  moved.  Typically,  heavy,  slow  ions  require  very  long  computational 
time  till  the  ion  distribution  has  converged.  To  speed  up  this  computation,  we  use 
hypothetically  light  ions,  which  have  the  same  mass  as  electrons.  The  use  of  such  a 
light  ion  is  justified  by  the  fact  that  we  are  only  interested  in  the  ion  density,  and  the 
assumption  that  no  ion  is  absorbed  by  the  tether  and  hence  that  the  ion’s  distribution 
function  is  Maxwellian.  Prom  this  assumption,  we  can  show  that  the  ion  density  does 
not  depend  on  the  ion  mass.  That  is, 


(3.44) 


This  fact  enables  us  to  use  the  very  light  ion,  with  which  the  computation  gives  the 
ion  density  in  the  computational  domain.  We  should  note  here  that  the  trajectory 
of  light  ions  is  different  from  that  of  real  ions,  but  the  corresponding  density  is  the 
same.  Actually  with  these  light  ions  (ra*  =  me ),  we  can  get  the  same  result  as  with 
other  ion  mass,  say  ra*  =  1840rae,  and  that  more  than  5  times  as  fast. 


3.4  Results 

In  this  section,  we  show  the  results  from  our  simulation.  First,  to  check  our  simulation 
qualitatively,  we  show  and  examine  the  field  quantities.  In  Figures  3-8,  3-9,  and  3-10, 
instantaneous  charge  densities  for  the  cases  of  the  Debye  ratio  equal  to  1,  2  and  5 
are  shown.  As  we  would  expect,  in  the  vicinity  of  the  tether,  there  is  a  region  where 
electrons  are  dominant  and  quasi-neutrality  no  longer  prevails.  This  region  is  the 
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“sheath” . 


From  this  electric  charge  density,  we  solve  by  Successive  Line  OverRelaxaion 
(SLOR)  to  obtain  the  electric  potential  distributions,  which  are  shown  in  Figure 
3-11,  3-12  and  3-13  for  the  same  cases.  These  figures  clearly  illustrate  the  sheath 
region  in  the  vicinity  of  the  tether.  The  electric  potential  there  is  positive,  indicating 
the  non-neutrality. 

Next  we  examine  the  quantitative  result  from  the  simulation.  Figure  3-14  shows 
the  history  of  current  collected  by  the  tether  and  the  corresponding  value  from  the 
results  of  Laframboise.  After  some  perturbation,  the  observed  current  oscillates  just 
below  the  Laframboise ’s  value.  This  oscillation  is  attributed  to  the  small  timestep  and 
the  small  number  of  particles  used  in  this  simulation.  The  reason  why  the  current 
collection  in  the  case  of  =  1  has  a  larger  amplitude  than  the  others  is  that  the 
surface  area  of  the  tether  is  smaller  than  that  of  the  others  and  thus  one  particle 
difference  becomes  more  significant  to  the  current  density  calculation. 

The  consistent  negative  bias  in  the  result  (about  7%)  is  probably  due  to  the  grid 
distortion,  as  noted  before.  But  more  work  needs  to  be  done  to  verify  this. 

Finally,  the  figure  3-15  shows  the  comparison  of  the  current  collection  obtained 
in  the  simulation  with  the  analytical  values  for  different  cases  of  the  Debye  ratio. 
From  this  figure,  the  PIC  method  is  verified  to  give  the  same  trend  of  the  analytical 
solution.  In  the  figure  3-15,  current  collected  is  normalized  by  the  random  thermal 
current  given  by  _ 


1  = 


Srcoolgl, 


2tt  m^e 
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Figure  3-15:  Current  collection  vs.  l/£p.  OML  current  (dashed),  the  exact  value 
(solid)  and  computed  value  (error  bars) 
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Chapter  4 


Conclusion  &;  Further  Work 


In  this  chapter,  we  conclude  the  possibility  and  perspective  of  PIC  method  to  be 
applied  to  the  space  engineering,  based  on  what  we  discussed  in  previous  chapters. 
First,  we  re-state  the  result  from  the  PIC  simulation. 

(1)  PIC  method  qualitatively  shows  the  same  trend  of  the  current  collec¬ 
tion  with  the  analytical  solution 

Figure  3-14  illustrates  the  exact  current  collection  as  a  function  of  the  Debye  ratio  and 
the  value  obtained  from  our  simulation.  As  the  Debye  ratio  decreases,  the  amount  of 
current  collected  by  the  tether  gets  closer  to  the  limiting  value.  This  limiting  value 
is  the  one  shown  by  equation  (2.13).  And  as  the  ratio  increase,  the  current  collection 
decreases.  This  is  because  of  the  bump  in  the  effective  potential  which  prevents  the 
particles  from  reaching  the  surface  of  the  tether  even  if  the  particle  has  an  enough 
energy  to  reach  the  surface. 

(2)  PIC  method  quantitatively  shows  the  current  collection  with  7%  error 
from  the  analytical  value  for  all  cases  of  Debye  ratio 

As  illustrated  in  Figure  3-14,  the  amount  of  current  collected  by  the  tether  oscillates 
just  below  the  analytical  value.  The  cause  of  the  error  can  be  attributed  to  the 
distortion  of  the  grid,  as  noted  before.  To  avoid  this  problem,  we  plan  to  use  different 
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methods  for  the  interpolation.  One  possibility  is  to  use  the  triangular  grid  with 
the  area-weighting  function  analogous  to  one  we  use  in  this  thesis.  The  use  of  the 
triangular  grid  does  not  need  the  linear  interpolation,  therefore  we  can  expect  to  get 
rid  of  the  error  caused  by  the  distortion  of  the  grid. 

(3)  The  current  collection  code  requires  non-zero  potential  as  an  outside 
boundary  condition. 

The  partial  absorption  of  electrons  requires  non-zero  potential  region,  “pre-sheath” 
outside  the  sheath.  The  pre-sheath  extends  very  far.  To  run  the  simulation  on  a 
relatively  small  grid,  we  need  to  know  the  local  potential  at  the  computational  outside 
boundary.  This  potential  was  calculated  from  the  equation  of  quasi-neutrality.  The 
quasi-neutrality  equation  relates  the  local  potential  to  the  instantaneous  outbound 
electron  density.  We  calculated  the  density  from  the  outbound  electron  flux. 

(4)  The  current  collection  code  can  reduce  the  computational  time  by 
using  a  hypothetically  light  ion  mass. 

Since  we  are  only  interested  in  the  ion  density,  not  in  the  ion  trajectories,  we  can 
use  a  very  light  ion  mass.  The  reasonable  assumption  that  no  ion  is  absorbed  by  the 
tether  allows  us  to  use  a  very  small  ion  mass.  Since  the  ion  density  does  not  depend 
on  the  ion  mass,  when  the  ions  have  Maxwellian  distribution,  the  light  ion  becomes 
a  critical  factor  to  speed  up  the  computation  considerably. 

Based  on  this  PIC  method,  we  plan  to  consider  (1)  the  case  of  flowing  unmagne¬ 
tized  plasma,  and  (2)  the  case  of  magnetized  plasma  as  seen  in  the  space  engineering 
application.  For  either  case,  the  tether  potential  is  still  high  and,  therefore,  we  can 
still  use  the  light  ion  mass.  As  for  the  symmetry  of  the  phenomena,  we  have  not  taken 
into  account  so  far  any  symmetric  conditions  except  for  the  outside  boundary  poten¬ 
tial,  which  is  taken  as  an  average.  This  only  symmetric  condition  can  be  removed  by 
applying  local  potential  values.  Therefore,  practically  this  method  does  not  need  any 
symmetry  conditions.  For  future  applications,  where  the  symmetry  condition  is  no 
longer  valid,  we  can  apply  our  PIC  method  without  a  substantial  modification. 
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